REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions, .searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  °r  °,hel[  aspect  ofthis 

collection  of  information, including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  DJri“°^!fnX"i!,0rrwo704  0188)  Washincrton  DC  20503^  ” 

Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  Dl  2UbUJ. _ 

1.  AGENCY  USE  ONLY  (Leave  blank!  I  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

25.Oct.99  THESIS 

4.  TITLE  AND  SUBTITLE - -  5.  FUNDING  NUMBERS 

MULTIRESOLUTION  IMAGE  FUSION  OF  THEMATIC  MAPPER  IMAGERY 
WITH  SYNTHETIC  APERTURE  RADAR  IMAGERY 


6.  AUTHOR(S) 

2D  LT  MEEK  THEODORE  R 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

UTAH  STATE  UNIVERSITY 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

THE  DEPARTMENT  OF  THE  AIR  FORCE 
AFIT/CIA,  BLDG  125 
2950  P  STREET 
WPAFB  OH  45433 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

FY99-339 


11.  SUPPLEMENTARY  NOTES 

12a.  DISTRIBUTION  AVAILABILITY  STATEMENT 

Unlimited  distribution 

In  Accordance  With  AFI  35-205/AFIT  Sup  1 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  ( Maximum  200  words) 

1999 

Distribution  Unlimited 

1117  076  i 

14.  SUBJECT  TERM 


15.  NUMBER  OF  PAGES 
114 


16.  PRICE  CODE 


17  SECURITY  CLASSIFICATION  I  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT  ABSTRACT 


DTIC  QUALITY  INSPECTED  4 


Standard  Form  298  (Rev.  2-89)  (EG) 

Prescribed  by  ANSI  Std.  239.18 

Designed  using  Perform  Pro,  WHS/DIOR,  Oct  94 


MULTIRESOLUTION  IMAGE  FUSION  OF  THEMATIC  MAPPER  IMAGERY  WITH 
SYNTHETIC  APERTURE  RADAR  IMAGERY 

by 

Theodore  R.  Meek 

A  thesis  submitted  in  partial  fulfillment 
of  the  requirements  for  the  degree 

of 

MASTER  OF  SCIENCE 
in 

Electrical  Engineering 


UTAH  STATE  UNIVERSITY 
Logan,  Utah 


Copyright  Theodore  R.  Meek  1999 
All  Rights  Reserved 


Ill 


ABSTRACT 

Multiresolution  Image  Fusion  of  Thematic  Mapper  Imagery 
with  Synthetic  Aperture  Radar  Imagery 

by 

Theodore  R.  Meek,  Master  of  Science 
Utah  State  University,  1999 

Major  Professor:  Dr.  Doran  J.  Baker 
Department:  Electrical  and  Computer  Engineering 

This  study  was  designed  to  demonstrate  the  feasibility  of  applying  multiresolution 
image  fusion  techniques  to  synthetic  aperture  radar  (SAR)  and  Landsat  imagery.  This 
was  accomplished  through  the  development  and  application  of  image  fusion  software  to 
test  images,  for  the  purpose  of  comparing  results,  to  show  that  information  from  more 
than  three  bands  can  be  used  to  study  surface  and  subsurface  features  in  a  single  image. 
The  test  images  were  fused  using  six  image  fusion  techniques  that  are  the  combinations 
from  three  types  of  image  decomposition  algorithms  (ratio  of  low  pass  [RoLP]  pyramids, 
gradient  pyramids,  and  morphological  pyramids)  and  two  types  of  fusion  algorithms 
(selection  and  hybrid  selection  and  averaging).  Based  upon  the  composite  images  formed 
by  fusing  the  test  images,  this  study  concludes  that:  small  details  in  city  areas  make 
morphological  pyramids  ineffective,  selection  forms  of  fusion  do  not  effectively  combine 
the  data,  RoLP  and  gradient  pyramids  with  hybrid  fusion  produce  the  best  results,  and 
optimum  pyramid  depth  is  dependent  upon  the  size  of  detail  in  the  images.  (122  pages) 
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CHAPTER  I 


INTRODUCTION 

A.  Remote  Sensing  and  Earth  Images 

Remote  sensing  is  the  technique  of  using  sensors  located  a  considerable  distance 
from  the  information  source  to  obtain,  process,  and  interpret  data  to  obtain  information 
and  perspective  about  the  source  itself.  Terrestrial  remote  sensing  has  become  popular  in 
the  last  two  decades  because  it  makes  possible  easy  mapping,  monitoring,  and  study  of 
changes  on  or  even  beneath  the  Earth’s  surface.  Images  obtained  from  remote  sensors 
have  been  used  for  the  study  of  ice  sheets  [1,2],  locating  villages  in  third  world  countries 
[3],  locating  gold,  vegetation  mapping,  land  use  studies,  geological  mapping,  watershed 
studies,  urban  planning,  and  much  more  [4,  5]. 

A  widely  used  source  of  remotely  sensed  data  is  the  Landsat  satellite  series. 
Landsat  uses  the  optical  and  near  infrared  (IR)  wavelengths.  Low  spatial  resolution  (83 
meters)  multispectral  scanner  instrument  images  from  the  first  Landsat  satellite,  Landsat- 
1,  have  been  available  since  1972  [6].  Although  Landsat  images  are  an  excellent  source 
of  information,  the  images  can  be  quite  expensive.  High  spatial  resolution  (30  meters) 
images  from  Landsat-5’s  Thematic  Mapper  (TM)  instrument  can  be  obtained  through 
Space  Imaging  EOS  AT  but  cost  thousands  of  dollars  per  scene.  A  less  expensive  source 
is  the  EROS  Data  Center  that  sells  images  over  10  years  old  for  $200  to  $400  per  scene 
[6].  There  have  been  five  Landsat  satellites  successfully  launched  by  the  United  States; 
Landsat-6  was  lost,  and  Landsat-7  is  scheduled  for  launch  in  April  1999. 
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Another  source  of  remotely  sensed  data  is  synthetic  aperture  radar  (SAR).  SAR 
images  of  the  Earth’s  surface  have  been  available  since  the  late  1970s.  The  first  U.S. 
SAR  satellite  was  Seasat,  which  was  launched  in  1978  to  collect  data  on  sea-surface 
winds,  sea  ice  features,  and  ocean  topography  [4].  The  Space  Shuttle  has  flown  three 
SAR  instruments:  SIR-A  in  1982,  SIR-B  in  1985,  and  SIR-C  in  1994.  SAR  has  become 
increasingly  popular  in  the  1990s,  most  likely  because  of  new  image  processing 
techniques  that  incorporate  accurate  correction  of  acquisition  errors  such  as 
foreshortening.  Several  countries  have  launched  SAR  satellites;  these  include  the 
European  Space  Agency’s  ERS-1  and  ERS-2  launched  in  1991  and  1996,  Japan’s  JERS-1 
launched  in  1992,  and  Canada’s  Radarsat  launched  in  1995  [4]. 

SAR  has  a  distinct  advantage  over  passive  sensors,  such  as  the  ones  in  Landsat 
satellites,  in  that  SAR  does  not  rely  on  solar  illumination.  In  fact,  SAR  can  gather  data  in 
both  sunlit  and  dark  regions  as  well  as  in  regions  covered  by  clouds,  fog,  or  pollution. 
SAR  can  even  penetrate  several  near-surface  materials  such  as  sand,  foliage,  and  forest 
canopies  [4,  7].  This  capability  makes  SAR  an  excellent  companion  for  study  with 
optical  sensors. 

B.  Viewing  Techniques  and  Challenges 

Most  remotely  sensed  data  come  in  bands  corresponding  to  the  wavelength  used 

to  acquire  the  images.  Typical  Landsat  TM  images  have  seven  bands,  and  SAR  images 
can  have  several  bands.  In  order  to  view  a  single  band,  there  are  two  main  options, 
representing  the  image  as  a  grayscale  image  or  as  a  pseudocolored  image. 
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To  view  a  single  band  using  the  grayscale  option,  each  sensor  sample,  or  pixel, 
will  be  assigned  a  value,  usually  between  0  and  255.  If  the  sensor  does  not  receive  or 
store  the  samples  on  a  0  to  255  scale,  the  conversion  is  simple:  The  highest  value  in  the 
range  corresponds  to  255;  the  lowest  corresponds  to  0,  and  all  other  values  are 
interpolated  and  quantized  to  fit  the  0  to  255  scale.  Color  pictures  are  made  up  of  three 
grayscale  images,  a  red  image,  a  green  image,  and  a  blue  image  (RGB),  commonly  called 
RGB  values.  For  example,  the  expression  RGB  =  (1,  4, 2)  would  mean  that  the  red  pixels 
are  band  1,  green  pixels  are  band  4,  and  blue  pixels  are  band  2.  In  a  grayscale  image,  it  is 
often  difficult  to  distinguish  between  pixels  of  similar  value  because  they  seem  to  blend 
together.  Pseudocoloring  is  used  because  the  human  eye  readily  detects  contrast  in 
images  [8,  9].  The  pseudocoloring  technique  is  to  assign  each  pixel  intensity  a  red, 
green,  and  blue  corresponding  color  level.  This  scale  helps  viewers  see  large-scale 
changes  across  the  image.  The  color  level  assigned  to  each  pixel  value  can  be  arbitrary; 
however,  the  following  is  a  common  scale:  Pixel  intensity  from  0  to  47  changes  from 
black  to  violet;  intensity  from  48  to  79  changes  from  violet  to  blue;  and  additional 
increments  of  47  go  from  blue  to  cyan  to  green  to  yellow  to  red  to  white  [5].  White  is 
still  the  most  intense  pixel,  and  black  is  the  least  intense. 

To  view  multiple  bands,  several  techniques  have  been  implemented.  Nonetheless, 
each  technique  is  based  upon  the  same  principle:  Assign  one  band  to  each  pixel  color 
red,  green,  and  blue.  This  limits  the  viewing  of  an  image  to  three  bands;  however,  it 
allows  full  color  images  of  the  data.  Typical  Landsat  images  are  viewed  with  color 
combinations  that  show  the  bands  of  interest  to  the  study;  different  bands  show  different 
properties  on  the  Earth’s  surface.  More  details  on  what  specific  Landsat  bands  show  will 
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be  given  in  Chapter  II.  Other  than  a  simple  band  being  used  to  represent  the  color,  ratios 
of  bands  have  also  been  used.  For  example,  red  could  be  band  7  divided  by  band  4;  green 
could  be  band  2  divided  by  band  1,  and  blue  could  be  band  7  divided  by  band  3  [10]. 

SAR  images  with  more  than  one  band  can  also  be  viewed  this  way;  for  example,  red 
could  be  band  L,  and  green  and  blue  could  both  be  band  C  [5]. 

Composite  images  can  be  formed  if  proper  image  registration  has  been  achieved. 
This  usually  requires  extensive  processing  of  the  SAR  images  to  correct  for 
foreshortening.  This  is  done  by  assigning  bands  from  both  sources  to  the  red,  green,  and 
blue  pixel  values.  An  example  would  be  to  use  two  Landsat  bands  for  red  and  green  and 
a  SAR  band  for  blue.  This  technique  has  proven  extremely  useful  in  determining  land 
usage;  however,  this  technique  is  still  limited  because  only  three  bands  can  be  viewed 
simultaneously. 

C.  Multiresolution  Image  Fusion 

Combining  two  images  to  form  a  composite  image  is  not  unique  to  remote 
sensing.  Several  applications  have  use  for  image  fusion,  i.e.,  the  merging  of  two  images 
to  create  one.  Image  fusion  can  be  used  for  computer  vision,  aiding  pilots  to  land  in 
inclement  weather,  compressing  source  image  data,  medical  diagnosis  through  fusing  CT 
and  MR  images,  ease  of  viewing  to  reduce  human  operator  work  load,  etc.  [11,  12,  13] 

The  simplest  approach  to  image  fusion  is  to  set  the  composite  pixel  value  equal  to 
the  average  of  the  source  values.  This  is  not  the  optimal  approach,  however,  because 
features  that  appear  in  one  image  and  not  the  other  are  shown  in  the  composite  at  a 
reduced  contrast  or  are  superimposed  on  features  from  the  other  image. 
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The  key  to  successful  image  fusion  is  to  create  a  composite  image  that  retains  all 


useful  information  contained  in  the  source  images  without  introducing  artifacts  that  can 
interfere  with  analysis  and  interpretation  [11,  12].  Figure  1  represents  the  approach  taken 
in  current  literature,  a  multiresolution  pyramid  [1 1-17].  In  this  approach,  the  source 
images  are  filtered  and  reduced  repetitively,  forming  a  pyramid.  In  the  figure,  Ao  and  Bo 
are  the  original  source  images  to  be  fused;  Ai  and  Bi  are  reduced- 


Composite  Pyramid 


Fig.  1.  Data  flow  diagram  for  pyramid  fusion  process 
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filtered  versions  of  Ao  and  Bo.  Likewise,  A2  and  B2  are  reduced-filtered  images  of  Ai 
and  Bi.  This  process  will  be  explained  in  more  detail  in  Chapter  HI.  Each  level  of  the 
pyramid  is  then  fused,  using  a  fusion  algorithm  to  create  a  composite  pyramid.  The 
fusion  process  can  be  either  complex  or  simple.  (Fusion  processes  are  discussed  in 
Chapter  IV.)  A  composite  image  is  generated  from  the  composite  pyramid,  which  will 
include  all  important  information  from  the  source  images  and  no  new  artifacts.  Figure  1 
is  a  data  flow  diagram  of  the  multiresolution  fusion  process.  The  square  boxes  represent 
the  levels  of  the  pyramid,  or  images.  The  red  circles  represent  functions  that  fuse  the  two 
incoming  images  into  one  composite  outgoing  image.  Although  no  function  appears 
between  levels  of  the  pyramids,  it  should  be  understood  that,  when  constructing  a 
pyramid  from  the  source  image,  there  is  a  function  that  filters  and  subsamples  the  current 
image  to  generate  the  higher  level.  The  green  boxes  represent  the  composite  pyramid; 
however,  Co  is  not  the  final  composite  image.  The  final  image  is  obtained  by 
reconstructing  the  composite  pyramid. 

D.  Study  Approach 

Analysis  of  Landsat  and  SAR  images  is  conventionally  accomplished  by 
assigning  particular  bands  or  ratios  of  bands  to  red,  green,  and  blue  pixels.  By  choosing 
specific  bands  or  ratios  of  bands,  analysis  of  particular  properties  can  be  made  easier.  In 
several  cases  it  is  possible  to  generate  images  from  different  bands  that  show  the  same 
properties.  A  popular  way  to  analyze  Landsat  and  SAR  images  of  the  same  areas 
concurrently  is  done  by  overlaying  bands  from  one  source  with  bands  of  another  source 
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after  the  images  have  been  properly  registered  [1,  2,  3,  5,  7],  for  example,  having  red  and 
green  pixels  represent  two  Landsat  bands  and  blue  pixels  represent  a  SAR  band. 

Current  image  fusion  techniques  can  create  a  composite  image  that  contains  all  of 
the  important  information  in  the  source  images  without  adding  new  artifacts  to  interfere 
with  image  analysis.  Since  a  human  observer  cannot  reliably  integrate  visual  information 
by  viewing  multiple  images  separately  and  consecutively  [17],  image  fusion  makes  it 
possible  for  the  observer  to  view  one  image  that  contains  all  of  the  information  contained 
in  the  source  images.  The  present  study  includes  programming  the  image  fusion 
techniques  described  in  current  literature  and  applying  them  to  Landsat  and  SAR  images 
to  diagnose  the  usability  of  each  fusion  technique  as  applied  to  remotely  sensed  imagery. 
The  results  and  composite  images  are  given  in  Chapter  V. 

E.  Study  Objectives 

The  first  objective  of  this  study  was  to  compare  current  image  fusion  techniques 
and  diagnose  their  effectiveness  in  fusing  SAR  and  Landsat  imagery  based  upon  the 
composite  image  results  obtained  through  the  fusion  of  real  SAR  and  Landsat  images  and 
the  time  required  to  generate  the  composite  image.  This  was  accomplished  by  applying 
the  fusion  techniques  to  several  test  images,  recording  the  time  required  for  each 
technique  to  produce  the  composite  image,  and  visually  inspecting  the  image  to  observe 
composite  image  quality. 

The  second  objective  was  to  combine  test  images  by  applying  the  fusion 
techniques  to  real  SAR  and  Landsat  image  data  to  obtain  useful  information  from  more 
than  three  remote  sensor  bands.  This  was  accomplished  by  creating  test  images  designed 
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to  show  specific  surface  features  and  then  fusing  the  test  images  to  verify  that  the  surface 
features  have  been  successfully  fused. 

The  third  objective  of  this  study  was  to  determine  the  optimum  pyramid  depth  for 
the  fusion  of  remotely  sensed  imagery.  This  was  determined  by  examining  the  same  set 
of  test  images  fused  over  a  range  of  pyramid  depths. 
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CHAPTER  II 

LANDSAT  AND  SAR  IMAGES 

A.  Landsat  Image  Characteristics 

The  Landsat  images  used  in  this  study  were  acquired  by  Landsat-4’s  TM  sensor 
recorded  on  25  September  1994  in  path  number  38,  row  number  31.  The  images  were 
purchased/requested  by  the  Rocky  Mountain  NASA  Space  Grant  Consortium 
(RMNSGC)  and  have  been  used  in  previous  studies  by  both  David  Oliver  in  1996  [5]  and 
Bill  Pfaffin  1997  [18]. 

Landsat  TM  images  do  not  have  significant  geometric  distortions  because  the 
viewing  angle  of  the  sensors  is  nearly  perpendicular  to  the  Earth’s  surface.  In  TM 
images,  each  pixel  is  stored  as  an  8-bit  number  and  represents  either  a  30x30  meter  or 
120x120  meter  terrestrial  area.  The  Landsat  satellite  is  in  low  Earth  orbit  with  an  altitude 
of  705  km  and  a  16  day,  233  orbit  cycle.  The  TM  sensor  has  seven  bands;  bandwidths, 
resolutions,  and  percents  of  transmission  for  the  TM  bands  are  listed  in  Table  I.  The 
percent  of  transmission  for  each  band  is  the  minimum  value  of  transmission  within  the 
corresponding  bandwidths.  Values  of  transmission  were  taken  from  [19]. 

Each  of  the  seven  bands  was  designed  to  maximize  detecting  and  monitoring  of 
different  types  of  Earth  resources  [20].  Band  1  is  in  the  visible  spectrum  corresponding 
to  blue  light;  it  penetrates  water  for  bathymetric  mapping  along  coastal  areas  and  is  used 
for  soil/vegetation  differentiation  and  for  distinguishing  forest  types.  Band  2  is  also  in 
the  visible  spectral  region  corresponding  to  green  light;  it  is  used  to  detect  green 
reflectance  from  healthy  vegetation.  Band  3  is  in  the  visible  spectrum  corresponding  to 
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TABLE  I 

LANDSAT  TM  BAND  INFORMATION 


Band 

Bandwidth  (|im) 

Resolution  (m) 

%  Transmission 

1 

0.45-0.52 

30 

100 

2 

0.53-0.61 

30 

95 

3 

0.62-0.69 

30 

100 

4 

0.78-0.91 

30 

84 

5 

1.57-1.78 

30 

98 

6 

10.42-11.66 

120 

88 

7 

2.08-2.35 

30 

98 

red  light;  it  is  designed  to  detect  chlorophyll  absorption  in  vegetation.  Band  4  is  in  the 
near-infrared  and  is  used  for  detecting  reflectance  peaks  in  healthy  green  vegetation  and 
for  detecting  water-land  interfaces.  Bands  5  and  7  respond  in  the  near-infrared  and  are 
used  for  vegetation  and  soil  moisture  observations  as  well  as  for  discriminating  between 
rock  and  mineral  types.  Band  6  is  selected  for  the  atmospheric  spectral  window  and  is  in 
the  thermal- infrared  spectrum,  designed  to  assist  in  thermal  mapping  and  in  soil  moisture 
and  vegetation  studies  [5, 20]. 

B.  SAR  Image  Characteristics 

Synthetic  aperture  radar  is  based  on  technology  that  uses  the  motion  of  the  sensor 
to  synthesize  an  antenna  aperture  larger  than  the  physical  antenna  [4].  This  technique 
provides  an  enhanced  spatial  resolution  imaging  capability  that  is,  to  the  first  order, 
independent  of  sensor  altitude.  SAR  functions  by  sending  a  microwave  signal  at  the 
target  and  recording  the  returned  signal.  Being  an  active  rather  than  a  passive  system, 
SAR  has  the  advantage  of  being  independent  of  solar  illumination  of  the  Earth’s  surface. 
SAR  return  signals  are  not  inhibited  when  the  Earth’s  surface  is  covered  by  clouds,  fog. 
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haze,  or  smoke;  in  fact,  SAR  can  even  penetrate  the  Earth’s  surface  to  return  information 
about  subsurface  features  [7]  or  can  penetrate  forest  canopies  that  obscure  the  ground  [4]. 
For  these  reasons,  SAR  can  produce  radar  terrain  maps  that  show  the  presence  of  major 
subjacent  features  not  shown  in  Landsat  images  at  the  same  spatial  resolution  [7]. 

SAR  images  are  produced  from  the  signals  reflected  from  the  Earth’s  surface. 

The  strength  of  the  return  signal  is  dependent  upon  both  the  surface  characteristics  and 
the  radar  system.  Table  II  gives  a  list  of  parameters  that  affect  the  return  signal.  When 
the  surface  roughness  is  at  the  same  scale  as  the  radar  wavelength,  signal  interference 
occurs  resulting  in  a  speckled  image.  The  look  angle  can  also  cause  problems,  such  as 
geometric  distortion  and  the  appearance  of  near  range  pixels  being  brighter  than  far  range 
pixels,  in  the  image  [1]. 

There  are  numerous  possible  SAR  wavelength  bands;  however,  unlike  the 
Landsat  TM  imagery,  different  SAR  bands  are  identified  by  letters.  The  surface 
properties  proposed  for  study  should  be  used  to  determine  the  band  used;  different  bands 
show  different  properties.  When  viewing  single-band  SAR  images,  a  user  can  make  use 


TABLE  II 

PARAMETERS  AFFECTING  SAR  RETURN  SIGNALS 


Ground  Parameters 

Radar  System  Parameters 

Terrain  Slope  and  Texture 
Surface  Roughness 

Complex  Dielectric  Constant 
Terrain  Feature  Orientation 
Shadowing 

Frequency,  Wavelength 

Radar  Look  Angle/Depression  Angle 
Antenna  Look  Direction 

Polarization 

Signal-to-Noise  Ratio 
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of  the  same  grayscale  or  pseudocoloring  techniques  discussed  for  TM  images  in  Chapter 
I.  For  SAR  images,  the  highest  pixel  values  represent  areas  of  saturated  radar  return; 
these  areas  would  appear  white  in  the  coloring  schemes  mentioned. 

Images  from  the  SIR-C  SAR  aboard  STS-68  are  used  in  this  study.  SIR-C 
operated  at  two  frequencies,  the  C  band  at  5.3  GHz  and  the  L  band  at  1.24  GHz.  Each 
band  has  five  images  that  correspond  to  different  polarization  modes.  The  five  modes  are 
designated  as  HH,  HV,  VH,  VV,  and  total  power.  The  first  letter  represents  the 
transmitted  signal,  and  the  second  letter  represents  the  return  signal;  H  designates 
horizontal  polarization,  and  V  stands  for  vertical  polarization.  Therefore,  the  VH  image 
would  be  the  data  acquired  from  a  vertically  polarized  transmitter  and  a  horizontally 
polarized  receiver.  Total  power  images  are  the  sum  of  the  horizontal  and  vertical 
transmissions  and  the  horizontal  and  vertical  receptions.  In  this  study,  the  L  and  C  band 
total  power  images  are  used.  These  images,  like  the  Landsat  images,  were  acquired  by 
the  RMNSGC  [5].  The  SIR-C  sensor  look  angle  is  35  degrees  off  nadir.  (Nadir  is 
straight  down.)  Each  pixel  covers  a  12.5  meter  by  12.5  meter  area. 

C.  Combined  Image  Characteristics 

The  capability  of  SAR  to  operate  in  either  light  or  dark  all-weather  environments, 
coupled  with  the  fact  that  radar  terrain  maps  made  from  SAR  data  show  the  presence  of 
major  subjacent  features  that  are  not  detectable  in  Landsat  pictures  at  the  same  scale, 
implies  that  SAR  would  be  an  excellent  complement  to  Landsat  images.  However,  the 
look  angle  of  the  SAR  sensor  results  in  large  geometric  distortions  in  the  resulting  image 
display.  Due  to  these  large  geometric  distortions,  some  previously  integrated  SAR  and 
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Landsat  studies  concluded  that  SAR  added  little  value  because  of  scene  registration 
problems.  New  SAR  image  processing  techniques  have  been  used  to  overcome  this 
problem,  and  SAR  images  can  now  be  registered  to  Landsat  images  with  errors  reduced 
to  the  order  of  one  pixel  [5].  Integrated  Landsat  and  SAR  images  have  been  used  to  aid 
in  geologic  mapping  [10,  21],  the  study  of  ice  sheets  [1],  and  for  land  use  studies  [5]. 

One  article  suggested  that  neither  the  Landsat  TM  nor  the  SIR-B  SAR  imagery  alone 
could  distinguish  the  land  use  patterns  in  Sudan,  while  an  integrated  SAR  and  Landsat 
approach  was  successful  [3]. 

In  the  present  study  only  images  that  have  already  been  processed  and  registered 
are  used.  Correction  of  geometric  distortions  of  the  SAR  images  and  SAR  and  Landsat 
image  registration  was  performed  by  David  Oliver  as  part  of  his  thesis  in  1996  [5]. 

Figure  2  and  Figure  3  are  the  Landsat-4  and  SAR  images  of  Cache  Valley  used  in  this 
study.  Figure  2  is  composed  of  Landsat-4  TM  bands  1 ,  4,  and  2  corresponding  to  the  red, 
green,  and  blue  pixel  values,  respectively.  Figure  3  is  composed  of  the  SAR  L  and  C 
total  power  bands.  The  red  pixels  correspond  to  the  L  band,  and  the  green  and  blue  pixels 
correspond  to  the  C  band.  Note  that  in  the  upper  right-hand  comer  of  Figure  3,  the  SAR 
image  contains  no  data;  this  is  due  to  the  correction  for  geometric  distortion. 

Corresponding  subsections  of  the  images  in  Figure  1  and  in  Figure  2  are  fused, 
using  the  multiresolution  techniques  described  in  current  literature.  The  usefulness  of 
these  fusion  techniques  as  applied  to  remote  sensed  imagery  will  be  assessed  in  Chapter 
V.  Different  band  combinations  can  be  used  to  show  different  surface/subsurface 
properties;  in  several  cases,  there  are  more  than  three  bands  that  can  be  used  to  show  the 
same  features.  To  more  effectively  analyze  multiple  images,  image  fusion  can  be  used. 


14 


Fig.  2.  Landsat  bands  1,  4,  2. 
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CHAPTER  ffl 

MULTIRESOLUTION  PYRAMIDS 

A.  Image  Structure  and  Multiresolution 

Pyramids 

This  study  will  use  three  techniques,  described  in  this  chapter,  to  represent  an 
image  over  a  range  of  scales.  Other  techniques  are  also  presented  here  because  they  are 
necessary  for  understanding  the  three  techniques  used  in  the  study.  The  techniques  used 
to  represent  the  images  over  a  range  of  scales  are  called  pyramids;  the  three  of  interest  are 
RoLP,  gradient,  and  morphological  pyramids. 

A  complete  image  description  can  be  obtained  by  studying  an  image  structure 
over  a  range  of  scales.  When  we  zoom  in  on  an  image,  we  clearly  see  the  substructure; 
however,  we  lose  the  clarity  of  the  outlines.  On  the  other  hand,  when  we  zoom  out  to 
look  at  the  entire  picture,  the  scene  loses  detail.  It  logically  follows  that  relevant  details 
of  an  image  can  be  observed  only  within  a  certain  range  of  spatial  resolution.  If  we  focus 
on  the  small  details,  we  lose  focus  of  the  big  picture;  on  the  other  hand,  if  we  zoom  out  to 
see  the  whole  picture,  it  is  difficult  to  discern  the  small  details.  Figure  4  illustrates  this 
point.  When  zoomed  out,  the  canopy  and  intake  of  the  F-15  can  be  easily  discerned; 
however,  when  zoomed  in  on  the  area  outlined  by  the  white  box  in  the  next  image,  the 
exact  position  of  the  canopy  and  intake  edges  are  not  as  easily  seen.  For  this  reason,  it  is 
desirable  to  represent  an  image  over  a  range  of  scales,  depending  upon  the  structural 
content. 

A  series  of  images  progressively  smaller  in  structural  content  can  be  created  by 
repetitive  application  of  a  processing  operator  with  a  progressively  increasing  scale.  This 
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Fig.  4.  Demonstration  of  necessity  for  multiple  resolutions. 


operator  would  eliminate  details  smaller  than  a  certain  size.  This  operator  acts  like  a 
filter,  just  as  sifting  gravel  through  screens  with  different  wire  spacings  filters  gravel  into 
different  groups,  dependent  upon  particle  size.  Repetitive  application  of  this  operator 
separates  the  image  into  scenes  with  different  resolution  of  detail.  By  reducing  the 
sample  frequency  and  increasing  the  filter  size,  a  hierarchical  relation  is  generated. 
Reducing  the  sample  frequency  is  the  same  as  subsampling  the  image.  A  pyramid  is  a 
sequence  of  images  in  which  each  image  is  a  filtered  and  subsampled  copy  of  its 
predecessor  [11].  The  term  “multiresolution  pyramid”  comes  from  the  relationship  where 
successive  levels  in  a  pyramid  are  reduced  resolution  copies  of  the  input  image  [11-17]. 

The  function  that  generates  the  next  level  of  the  pyramid  could  be  called 
REDUCE  since  both  the  resolution  and  sample  density  are  decreased.  REDUCE  would 
both  filter  and  subsample  the  image.  To  create  a  pyramid  starting  with  the  source  image 
asP0, 
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P*  =  REDUCE(P*./)  for  k  =  1,2 ,...n ,  (1) 

where  n  is  the  number  of  levels  in  the  pyramid. 

Pyramid  reconstruction  to  recover  an  image  from  its  pyramid  will  need  an 
EXPAND  function  because  each  level  differs  in  sample  density.  EXPAND  is  defined  as 
follows: 

Pk-i  =  EXPAND  (P*)  for  k=n-l,n- 2,...0,  (2) 

where  n  is  the  number  of  levels  in  the  pyramid.  Specific  details  of  the  EXPAND  and 
REDUCE  operators  are  dependent  upon  which  types  of  pyramids  are  used.  Techniques 
used  to  generate  pyramids  can  be  classified  into  two  types:  (1)  linear  and  (2) 
morphological.  Sections  B,  C,  D,  and  E  describe  linear  techniques,  and  Section  F 
describes  morphological  techniques. 

B.  Gaussian  and  Laplacian  Pyramids 

As  given  in  section  A,  each  image  in  a  pyramid  is  a  low-pass  filtered  and 
subsampled  copy  of  the  previous  image.  This  is  beneficial  for  our  purposes  because  it 
defines  the  image  features  in  a  set  of  images  with  different  spatial  resolutions  by  filtering. 
The  most  common  linear  low-pass  filter  used  for  pyramid  generation  is  convolution  with 
a  Gaussian  kernel.  Pyramids  formed  using  this  technique  are  commonly  referred  to  as 
Gaussian  pyramids  [11,  12,  14,  16,  17]. 

To  create  a  Gaussian  pyramid  from  an  image  I,  assume  Gt  is  the  kth  level  of  the 
pyramid.  The  bottom  level  of  the  pyramid,  Go,  equals  I;  and 


G*  =  [w  *  G*  -  /]  1 2 


for  k  >  0 , 


(3) 
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where  [ _ ]j,2  denotes  subsampling  by  two  (The  resultant  image  has  .5  the  x  and  y 

dimensions  of  the  source  image.),  *  is  the  convolution  operator,  and  w  is  the  Gaussian 
kernel.  All  of  the  elements  of  w  must  sum  to  equal  1 ;  and,  for  simplicity,  w  is  defined  as 
separable,  where  w  =  w  *  w  .  Convolution  with  the  w  matrix,  as  defined  below,  yields  a 
Gaussian-similar  curve;  with  w  defined,  w  can  be  derived.  The  w  matrix  defined  in  [12] 
is 
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which  yields  a  w  matrix  defined  as 


1  2  1 
2  4  2 
1  2  1 


v16y 


[12]- 


(5) 


It  should  be  noted  that  this  value  for  w  was  used  because  of  its  Gaussian-like  properties; 
however,  w  can  be  defined  differently  as  long  as  it  meets  the  criteria  mentioned  above. 

Low-pass  filtering  and  subsampling  can  be  combined  into  a  single  function  fitting 
the  description  of  the  REDUCE  operator  defined  in  the  previous  section.  The  REDUCE 
operator  can  be  implemented  as  a  function  where  each  pixel  in  the  new  image  is  set  to  the 
weighted  average  of  the  corresponding  area  of  pixels  in  the  previous  image.  This  process 
is  much  simpler  than  the  previous  process  of  convolution  and  subsampling  (one-quarter 
the  number  of  operations),  yet  it  yields  the  same  result.  This  is  because  the  value 
returned  by  the  weighted  average  is  the  same  as  the  value  returned  by  convolution. 
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From  the  weighted  average  method,  the  kth  level  of  the  pyramid  is  defined  on  a  pixel  by 
pixel  basis  as 

2 

Pk(i,  j)=  ^  w(m,  n)Pk  -  i(2i  +  m, 2  j  +  n)  for  k  >  0 .  (6) 

m,n=-2 

The  weighting  function  w(m,n )  is  separable, 
w(m,n )  =  w'(m)w'(n) . 


The  function  w’  is  normalized, 

^  w  (m)  =  1 ,  (7) 

m=- 2 

and  is  symmetric, 

w'ii)  =  w'(-i)  fori  =  0,1,2.  (8) 

All  of  these  constraints  are  satisfied  when  w’  is  defined  as  follows  [11]: 

w(0)=a;  (9) 

w,(l)  =  w'(-l)  =  l/4;  (10) 

w'(2)  =  w>'(-2)  =  1/4- o/ 2.  (11) 


The  value  of  a  used  in  this  study  was  0.4  because  of  its  Gaussian-like  properties 
[11,  12,  14].  The  reader  should  note  that  w  and  w  are  matrices  used  in  convolution, 
whereas  w  and  w’  are  functions  used  in  the  weighted  average. 

Taking  the  difference  between  two  levels  in  a  Gaussian  pyramid  results  in  a  band¬ 
pass  filtered  image.  In  order  to  take  the  difference  between  the  two  levels,  however,  the 
level  of  lower  resolution  must  be  expanded  to  the  same  size  as  the  image  at  the  higher 
resolution.  If  this  is  done  for  every  level  of  a  Gaussian  pyramid,  a  Laplacian  pyramid  is 

created.  Let  L  be  the  kth  level  of  the  Laplacian  pyramid,  defined  as 


L*  =  G*  -  4w  *[G*  +  /]  T  2 , 


(12) 
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using  the  convolution  notation,  where  [. . -]t2  indicates  upsampling  by  two  [12].  When 
upsampling,  rows  and  columns  of  zeros  are  added  between  the  existing  rows  and 
columns.  Convolution  with  w  interpolates  values  for  the  columns  and  rows  of  zeros. 

The  Laplacian  pyramid  can  also  be  defined  using  the  weighted  average  method  as 

L k  =  G*  -  EXPAND!  G*  +  /)  [  1 1 , 14].  (13) 

For  linearly  filtered  pyramids,  EXPAND  can  be  defined  as 

Pk(j,y')  =  4-  w(m,n)Pk- 1  — 1,  (14) 

where  only  integer  coordinates  contribute  to  the  sum.  It  is  noted  that  in  the  Laplacian 

pyramid  L.  =  G» ,  namely,  the  top  level  of  the  Laplacian  pyramid,  is  the  same  image  as 
the  top  level  of  the  Gaussian  pyramid. 

The  pyramid  algorithm  reduces  the  filter  band  limit  by  an  octave  from  level  to 
level.  The  Laplacian  pyramids  are  the  equivalent  of  a  band-pass  filter  with  a  bandwidth 
equal  to  the  span  between  octaves.  Filtering  an  image  has  the  effect  of  blurring.  This  is 
because  sharp  edges  are  composed  of  high-frequency  components,  and  the  low-pass  filter 
attenuates  the  high  spatial  frequencies  in  the  image.  Images  higher  in  the  pyramid 
structure  have  fewer  high-frequency  components  than  images  in  the  lower  part  of  the 
pyramid.  Therefore,  when  expanded  to  the  size  of  an  image  lower  in  the  pyramid,  the 
image  higher  in  the  pyramid  appears  blurry. 

Reconstructing  an  image  from  its  Laplacian  pyramid  is  straightforward.  The 
original  image  is  recovered  by  expanding  the  highest  level  of  the  pyramid  and  adding  that 
to  the  next  highest  level  repeatedly,  until  no  levels  are  left.  This  is  mathematically 
expressed  as 
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p> 

II 

r« 

(15) 

and 

G*  =  L*  +  EXPAND(G*  +  /) 

for  k  -n  —  l,n  —  2,...0. 

(16) 

A 

This  yields  an  exact  reconstruction  of  the  source  image.  Go  =  Go. 


C.  Filter-Subtract-Decimate 
Laplacian  Pyramids 

A  level  of  the  filter-subtract-decimate  (FSD)  Laplacian  pyramid  is  defined  as  the 
difference  between  the  Gaussian  level  and  the  filtered  copy  of  the  Gaussian  level  prior  to 
subsampling  for  the  next  Gaussian  level.  For  this  reason,  the  REDUCE  function  as 
defined  will  not  work  with  FSD  Laplacian  pyramids  because  REDUCE  filters  and 
subsamples  at  the  same  time.  Let  L*  be  the  kth  level  of  the  FSD  Laplacian  pyramid, 
which  is  mathematically  defined  as 

L*  =  G*-w*G*=[l-w]*G*,  (17) 


where  1  is  a  matrix  of  the  same  dimensions  of  w  where  all  values  are  zero  except  for  the 
center  value,  which  is  one.  In  the  case  of  a  3-by-3  matrix, 


1  = 


0  0  0 
0  1  0 
0  0  0 


(18) 


Inspection  of  the  process  of  forming  the  image  in  the  Laplacian  pyramid,  L* , 
reveals  that  the  Gaussian  image,  G*  ,  is  sequentially  convolved  with  w,  subsampled, 
upsampled,  convolved  with  w  again,  and  then  subtracted  from  itself.  Refer  to  equations  3 
and  12.  However,  if  the  resampling  steps  are  skipped,  the  results  are  only  slightly 
different.  We  can,  therefore,  say  that 


L*  =  G k  -  4w  *  [[w  *  Gt]  i  2]  r  2 «  G*  -  w  *  w  *  Gt 
=  [1  -  w  *  w]  *  Gt 
=  [l  +  w]*[l-w]*Gt. 
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(19) 

Hence,  the  FSD  Laplacian  pyramid  can  be  converted  to  a  Laplacian  pyramid  through  the 
following  conversion: 

L*  «[l+w]*L*.  (20) 

Pyramid  reconstruction  to  obtain  the  original  image  follows  by  converting  the 
FSD  Laplacian  pyramid  into  a  Laplacian  pyramid  and  then  reconstructing  the  Laplacian 
pyramid  [12].  Note  that  this  reconstruction  is  not  an  exact  replica  of  the  original. 

D.  Ratio  of  Low-Pass  Pyramids 

The  ratio  of  low-pass  (RoLP)  pyramid  gets  its  name  from  the  relationship  that 
exists  between  successive  levels  of  the  pyramid.  RoLP  pyramids  are  very  similar  to 
Laplacian  pyramids;  instead  of  taking  the  difference  between  levels  of  a  Gaussian 
pyramid,  the  RoLP  pyramid  takes  the  ratio  between  levels  of  a  Gaussian  pyramid.  The 
RoLP  pyramid,  Rt ,  is  mathematically  defined  as 


Rt  =  Gk 

EXPAND  (Gt  + ;) 

for  k  =n— l,n  —  2,...0 

(21) 

and 

O 

II 

(22) 

Every  level  in  the  RoLP  pyramid  is  the  ratio  of  two  successive  levels  in  the  Gaussian 
pyramid. 

A 

Let  Go  be  the  image  reconstructed  from  the  RoLP  pyramid.  The  reconstruction 
process  is  the  inverse  of  the  construction  process. 


Go  —  Rn  , 


(23) 


(24) 


and  Gt  =  Rt  EXPAND(Gt  +  /)  for  k  =  n  - 1,  n  -  2,...0 . 

A 

The  reconstruction  process  for  RoLP  pyramids  is  exact,  or  in  other  words,  Go  =  Go . 
There  is  also  a  contrast-enhanced  RoLP  pyramid  described  in  [17]  that  was  not  used  in 
this  study. 


E.  Gradient  Pyramids 

The  term  “gradient  pyramid”  is  a  misnomer,  because  a  gradient  pyramid  is 
actually  a  collection  of  four  pyramids.  Let  D tm  represent  the  kth  level  and  mth  orientation 
gradient  pyramid  image  for  an  image  I.  Dim  is  obtained  from  convolving  Gt  with  a 
gradient  filter  dm ,  also  called  an  “oriented  second  derivative  filter.” 

Dt™  =  dm  *  [Gt  +  w  *  Gt] ,  (25) 

where 


d/  =  [l  -l]; 


and 


d2  = 

di  = 

d4  = 


0  -1 
1  0 
-1' 

1 

-1  0 
0  1 


1 

& 


1 

vr 


(26) 

(27) 

(28) 
(29) 


Reconstruction  of  an  image  from  a  gradient  pyramid  requires  Laplacian  and  FSD 
Laplacian  pyramids  as  intermediate  steps.  Each  gradient  pyramid  level  Dim  is  converted 

to  a  corresponding  second  derivative  pyramid  or  oriented  Laplacian  level  Ltm . 
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An  FSD  Laplacian  pyramid,  L * ,  is  then  formed  by  summing  the  oriented  Laplacian 
pyramids. 

Lt  =  ££*«.  (31) 

m= 1 

Reconstruction  is  then  completed  by  following  the  same  steps  required  to  reconstruct  an 
FSD  Laplacian  pyramid  [12]  as  described  in  Section  C  of  this  chapter.  Since  the  FSD 

Laplacian  pyramid  is  used,  composite  image  reconstruction  is  approximate,  or  Go-Go. 

F.  Mathematical  Morphology  Pyramids 

Mathematical  morphology  is  the  examination  of  the  geometric  structure  of  an 
image  by  probing  its  microstructure  with  certain  elementary  forms  called  structuring 
elements  [15].  Linear  filters  alter  the  object  intensities  and,  therefore,  estimate  the 
location  of  contours.  Morphological  filters,  on  the  other  hand,  remove  details  without 
adding  grayscale  bias  and,  therefore,  are  nicely  suited  to  shape  extraction.  This  section 
contains  a  brief  summary  of  mathematical  morphology  and  its  operations.  For  more 
information  on  any  of  the  topics  in  this  section,  the  reader  is  referred  to  [1 1,  15, 22,  23]. 

Mathematical  morphology  is  based  upon  set  operations.  When  used  with  images, 
the  first  set  is  considered  to  be  the  image;  the  second  set  is  the  structuring  element.  The 
sets  are  representative  of  shapes  that  are  manifest  in  the  image.  Selection  of  the 
structuring  element  is  an  important  part  of  morphological  operations.  The  smaller  the 
structuring  element,  the  smaller  the  details  that  are  filtered  or  modified.  When  generating 
a  pyramid  by  successively  increasing  the  size  of  the  structuring  element  or  by  reducing 
the  image  size,  successively  larger  image  details  are  filtered  out,  once  again  leading  to  the 
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multiple  resolutions.  Also,  when  reconstructing  the  image,  the  positional  accuracy 
cannot  be  better  than  the  radius  of  the  structuring  element. 

The  building  blocks  of  morphology  are  the  dilation  and  erosion  operators. 
Dilation  is  a  transform  that  combines  two  sets  using  vector  addition  of  set  elements.  If  A 
and  B  are  sets,  the  dilation  of  A  by  B  is  the  set  of  all  possible  vector  sums  of  pairs  of 
elements,  one  from  A  and  one  from  B.  Dilation  is  denoted  by  the  ®  symbol  and  is 
mathematically  defined  as 

A  ©  B  =  {c  g  EN  |c  =  a  +  b  for  some  a  e  A  and  b  e  B} ,  (32) 

where  A  and  B  are  the  sets  in  Euclidean  N-space  (EN)  with  elements  a  and  b, 
respectively.  The  complexity  of  image  dilation  and  erosion  is  similar  to  that  of 
convolution.  To  easily  conceptualize  dilation,  picture  a  straight  line;  after  dilation,  this 
line  becomes  wider. 

Erosion  is  the  morphological  dual  of  dilation.  The  erosion  transform  combines 
two  sets  using  the  vector  subtraction  of  set  elements.  If  A  and  B  are  sets  in  Euclidean  N- 
space,  erosion  of  A  by  B  is  the  set  of  all  elements  x  for  which  x  +  b  €  A  for  every  be  B. 
Erosion  is  denoted  by  the  0  symbol  and  is  mathematically  defined  as 

A0B  ={xe  EN|x  +  be  Aforeverybe  B},  (33) 

where  A  and  B  are  the  sets  in  Euclidean  N-space  with  elements  a  and  b,  respectively. 
Once  again,  for  easy  conceptualization,  picture  a  thick  line;  everywhere  the  line  is  thick 
enough  to  contain  the  structuring  element,  it  becomes  thinner.  If  this  line  is  dilated  and 
then  eroded,  it  would  appear  unchanged  because  everywhere  dilation  made  it  wider, 


erosion  made  it  thinner. 
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Conceptually,  an  erosion  shrinks  image  features,  and  a  dilation  expands  image 
features.  Figure  5  shows  the  erosion  and  dilation  of  a  box.  Note  how  the  box  walls  look 
thinner  after  dilation  and  thicker  after  erosion;  if  this  seems  the  opposite  of  the  expected 
result,  reconsider  the  expected  result  focusing  on  white  areas  as  the  boxes.  Knowing  that 
morphology  is  based  on  set  theory,  consider  the  color  black  the  empty  set.  The  black 
walls  are  areas  where  no  information  exists.  (Pixel  intensity  =  0  for  red,  green,  and  blue 
pixels.)  Now  consider  that  the  color  white  is  the  full  set;  this  would  be  areas  containing 
all  three  colors:  red,  green,  and  blue.  Where  blue  is  seen,  only  blue  is  present;  where  red 
is  seen,  only  red  is  present;  and  in  areas  where  a  pinkish-purple  color  is  seen  (only  in 
Figure  5  (b)),  both  red  and  blue  are  present. 

The  source  image  has  four  areas  of  concern:  the  black  box  around  the  outside 
where  no  pixels  have  values  greater  than  zero;  the  white  area  between  the  black  and  blue 
boxes  where  red,  green,  and  blue  all  have  pixel  values  greater  than  zero;  the  blue  area 


(a)  (b)  (c) 

Fig.  5.  Example  of  erosion  and  dilation  operators  using  a  10x10  brick  structuring 
element,  (a)  Source  image,  (b)  Source  image  after  dilation. 

(c)  Source  image  after  erosion. 
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where  only  blue  pixels  have  values  greater  than  zero;  and  the  red  area  where  only  the  red 
pixels  have  values  greater  than  zero. 

After  dilation,  the  white  walls  become  thicker,  making  the  black  walls  thinner. 
Likewise,  both  the  red  walls  and  the  blue  walls  expand,  causing  a  new  area  where  both 
red  and  blue  values  exist  simultaneously  without  green.  Although  the  blue  and  red  box 
appear  smaller  in  size,  they  have  actually  grown;  we  do  not  see  their  growth  because,  in 
the  areas  where  they  have  grown,  other  pixel  values  exist  so  we  either  see  white  or  else 
the  combination  of  blue  and  red,  which  is  the  pinkish-purple  color  in  Figure  5  (b). 

When  erosion  occurs,  the  white  walls  become  thinner,  making  the  black  walls 
bigger.  When  the  color  boxes  are  viewed,  it  is  easy  to  see,  by  looking  at  the  size  of  the 
red  box,  that  erosion  makes  the  boxes  smaller.  The  black  box  around  the  red  box  in 
Figure  5  (c)  is  an  area  where  both  the  red  box  and  the  blue  box  have  been  eroded  and  no 
image  data  are  left,  hence  the  empty  set,  or  black.  For  more  information  about  the 
erosion  and  dilation  of  a  set  or  the  properties  associated  with  erosion  and  dilation,  refer  to 
[22,  23]. 

The  steps  for  creating  a  morphological  pyramid  are  the  same  as  creating  a 
Gaussian  pyramid.  To  generate  the  next  pyramid  level,  the  current  image  is  filtered  and 
subsampled. 

Morphological  filters  are  sequences  of  morphological  operations  that  have  special 
properties  with  respect  to  the  shapes  in  the  image.  Morphological  filters  are  idempotent 
and  increasing.  Idempotent  means  that  successive  applications  of  the  filter  leave  the 
result  unchanged  after  the  first  time  the  filter  has  been  applied,  similar  to  the  application 
of  a  linear  band-pass  filter  to  a  signal.  Increasing  means  that  the  operations  maintain 
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inclusion  relationships  on  the  images  they  transform  [15];  or  if  set  A  is  a  subset  of  set  B, 

then  A  filtered  by  F  is  a  subset  of  B  filtered  by  F. 

The  simplest  morphological  filters  are  the  opening  and  closing  transformations. 

An  opening  is  an  erosion  followed  by  a  dilation  and  is  denoted  by  the  o  symbol.  A 

closing  is  a  dilation  followed  by  an  erosion  and  is  denoted  by  the  •  symbol. 

Mathematically,  opening  and  closing  are  defined  as  follows: 

A  o  f  =  (f0A)  ©  A ;  (34) 

A  •  f  =  (f  ®  A)0A ,  (35) 

where  A  is  the  image  and  f  is  the  structuring  element.  These  filters  are  also  low-pass 

because  they  attenuate  high-frequency  fluctuations  between  the  set  and  its  complement. 

Openings  and  closings  are  considered  dual  operators  because  what  one  does  to  the 

foreground,  the  other  does  to  the  background.  To  ensure  that  the  foreground  and 

background  of  an  image  are  treated  the  same,  openings  and  closings  usually  follow  one 

another.  Morphological  filters  are  most  often  a  combination  of  openings  and  closings  or 

closings  and  openings.  For  this  study,  we  will  use  F  to  represent  a  closing-opening  filter: 

F  =  (A  •  f )  o  f  .  (36) 

With  morphological  filters  defined,  we  can  describe  the  generation  of 

morphological  pyramids  in  a  manner  similar  to  the  generation  of  Gaussian  pyramids.  Let 

I  represent  the  original  image.  The  base  of  the  pyramid,  Mo ,  needs  to  be  a 

morphologically  filtered  copy  of  I;  the  filter  used  to  generate  Mo  will  determine  what 

filters  are  to  be  used  for  pyramid  reconstruction.  We  will  again  define  the  REDUCE 

function,  which  will,  in  this  case,  morphologically  filter  the  image  and  then  subsample. 
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To  generate  the  next  level  in  the  pyramid,  REDUCE  will  be  applied  the  current  level.  Let 
Mo  be  the  base  of  the  pyramid,  and 

M  =  REDUCE(M  - /)  for  1  <  £  <  « ,  (37) 

where  n  is  the  depth  of  the  pyramid  and  REDUCE  is  F(M  - 1)  i  2.  From  here,  a 
difference  pyramid  similar  to  the  linear  Laplacian  pyramid  can  be  constructed  where 
D„  =M  (38) 

and 

D*  =  M  - EXPAND(M  +  /)  for  k=n~ln~ 2,...0 .  (39) 

Here,  EXPAND  is  defined  as  upsampling  followed  by  a  closing.  REDUCE  and 
EXPAND  can  use  any  filter;  however,  the  user  must  ensure  that  the  filters  used  in  the 
EXPAND  function  complement  the  filters  used  in  the  REDUCE  function  [22,  23].  If  the 
initial  filtering  of  the  source  image  I  is  a  closing,  then  a  minimal  reconstruction  should  be 
performed;  and  the  EXPAND  function  should  use  the  closing  filter.  On  the  other  hand,  if 
the  initial  filtering  of  the  source  image  I  is  an  opening,  then  a  maximal  reconstruction 
should  be  performed;  and  the  EXPAND  function  should  use  the  dilation  filter. 

Pyramid  reconstruction  is  once  again  straightforward.  The  reconstructed  image, 

A 

Mo ,  can  be  obtained  by  setting 

M,  =  D„  (40) 

and 

Mt  =  Dt  +  EXPAND(M  +  /)  for  k  =  n  - 1  ,n  -  2,...0 .  (41) 

A 

Mi  is  an  exact  reconstruction  of  the  pyramid  if  no  other  pyramid  processing  has  occurred 
to  either  the  morphological  pyramid  or  the  difference  pyramid. 
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CHAPTER  IV 

PYRAMID  FUSION  TECHNIQUES 
A.  Pyramid  Fusion  Overview 

This  study  will  use  three  of  the  pyramid  techniques  described  in  the  previous 
chapter  to  combine  or  fuse  two  source  images  into  a  single  composite  image.  In  order  to 
do  this,  we  must  define  a  way  to  fuse  two  pyramids  into  a  single  pyramid.  This  chapter 
will  introduce  the  fusion  techniques  described  in  current  literature  that  were  used  in  this 
study. 

Pyramids  are  simply  a  convenient  way  to  represent  an  image  over  a  range  of 
spatial  resolutions.  By  combining  the  images  at  each  level  of  the  pyramid,  the  composite 
image,  formed  by  pyramid  reconstruction,  will  have  consistency  over  all  resolutions. 

When  fusing  two  pyramids,  each  of  the  levels  of  the  pyramids  is  fused  into  a 
composite  level,  resulting  in  a  composite  pyramid.  Refer  to  Figure  1.  Once  the 
composite  pyramid  is  formed,  the  fused  image  of  the  source  images  is  generated, 
employing  the  pyramid  reconstruction  techniques  associated  with  the  technique  used  to 
generate  the  source  pyramids.  For  example,  if  Laplacian  pyramids  A  and  B  were 
generated  from  two  source  images,  the  composite  image  resulting  from  the  fusion  of 
pyramids  A  and  B  would  be  reconstructed  from  the  composite  pyramid  C,  using  the 
techniques  described  in  Chapter  III  to  reconstruct  an  image  from  its  Laplacian  pyramid. 
Each  level  of  the  composite  pyramid  is  defined  as 

C*  =  FUSE(A*,Bt)  for  k  =  n,  n  - 1,  n  -  2,...0 ,  (42) 
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and  n  is  the  number  of  levels  in  the  pyramid.  FUSE  is  a  function  that  converts  the  two 
images  into  the  composite,  using  a  fusion  algorithm. 


B.  Selection  and  Averaging  Fusion  Techniques 

When  combining  two  images,  the  untrained,  intuitive  approach  might  be  to 
average  the  pixels  from  the  source  images  to  obtain  the  value  for  the  corresponding  pixel 
in  the  composite  image.  This  approach  is  undesirable  because  features  that  appear  in  one 
image  and  not  the  other  will  show  up  in  the  composite  at  a  reduced  contrast  or  will 
appear  superimposed  on  features  from  the  other  image,  much  like  when  camera  film  is 
double  exposed.  To  avoid  the  fusion  problems  by  averaging,  the  composite  image  can  be 
obtained  by  selecting  pixels  from  either  of  the  source  images. 

The  simplest  pixel  selection  technique  for  RoLP  pyramids  is  to  use  the  local  area 
contrast  to  determine  which  pixel  to  select.  The  contrast  of  a  given  pixel  is  defined  as  the 
ratio  of  the  difference  in  pixel  intensity  to  area  intensity  to  the  area  intensity,  or  as 
follows: 


Contrast^,  j)  =  MJlzMJl  =  MJl  _ , , 
Lb(i,J)  U(i,j) 


(43) 


where  L  is  the  luminance  at  (i,  j),  or  simply  the  pixel  intensity,  and  Lb  is  the  background 
luminance  for  that  area.  Note  that  the  ratio  of  L  to  Lb  is  the  value  in  the  RoLP  pyramid; 
hence,  to  get  the  contrast  for  a  certain  pixel  value,  simply  subtract  one  from  the  RoLP 
value. 


The  human  eye  detects  contrast  very  well.  When  taking  advantage  of  this  fact, 
pixels  of  maximum  contrast  are  selected  from  each  source  image  to  form  the  composite. 
The  logic  behind  this  method  is  that  by  selecting  details  of  maximum  contrast,  the  fused 
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image  will  provide  better  details  for  human  analysis.  The  contrast  version  of  the  FUSE 
function  is  implemented  on  a  pixel-by-pixel  basis,  as  follows: 


C  (i,y)  = 


I.B  (ij). 


when  (|A  (i,  j)  - 1|  >  |B  (i,  7)  -  l|l 
otherwise 


[11,16,17].  (44) 


However,  another  reference  implemented  the  contrast  fusion  in  a  slightly  different  way, 
as  follows: 


C  =  \ 


A(i,  j) , 


when  (|A(i,  y)|  >  |B(i,  y)| 
otherwise 


[13]. 


(45) 


It  is  interesting  that  it  was  implemented  this  way  because,  without  subtracting  one  from 


the  RoLP  pyramid  value,  it  is  not  based  on  the  area  contrast  anymore. 


C.  Hybrid  Averaging  and  Selection  Fusion 

The  problem  with  contrast  fusion  is  that  it  is  susceptible  to  noise.  Noisy  images 
are  typically  of  higher  contrast.  Using  a  contrast  image  fusion  technique  would, 
therefore,  result  in  a  composite  image  with  more  noise.  The  combination  of  both 
selection  and  averaging  to  generate  a  composite  image  was  proposed  by  [12].  For  this 
method  to  work,  a  metric  is  necessary  to  indicate  when  to  use  selection  and  when  to  use 
averaging;  this  metric  is  called  “match.”  When  the  match  measure  is  a  value  that 
corresponds  to  using  selection,  it  is  necessary  to  know  from  which  image  to  select.  This 
requires  a  saliency  measure.  When  the  two  images  are  distinctly  different,  the  composite 
image  should  select  the  most  salient  component.  However,  when  the  two  images  are 
similar,  the  composite  image  should  contain  the  mean  of  the  two  source  images.  This 
technique  makes  it  possible  to  reduce  noise  without  double  exposure  artifacts  [12]. 
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A  good  measure  of  pattern  salience  is  pixel  intensity.  To  obtain  the  salience  for  a 
given  sample,  the  individual  pixel  intensity  can  be  used.  Alternatively,  the  average  pixel 
intensity  of  a  small  area  near  the  sample  can  be  used.  The  mean  pixel  intensity  can  be 
mathematically  defined  as  a  weighted  average  over  an  area,  p,  using  S*(i,  j)  as  the 
saliency  measure  of  the  pixel  at  (i,j)  for  the  image  k,  as  follows: 


m 

S*(i,  j)  =  X  P(*'’  /)p‘(f + J  +  /) . 


(46) 


where  m  can  range  from  0,  which  would  be  the  individual  sample,  to  2,  which  would 
include  the  5x5  area  surrounding  the  sample.  The  value  of  p(i',  /)  is  weight  of  the 
sample  within  the  area  or  neighborhood  p;  samples  closer  to  the  (i,  j )  position  have  higher 
weight  values.  The  function  p  serves  the  same  purpose  here  as  does  the  function  w’  in 
Chapter  III,  Section  B.  The  value  returned  by  P *(/  + 1,  j  +  f)  returns  the  pixel  value 
inside  the  image  at  location  (/,  j),  offset  by  i’  and  j’;  it  can  simply  thought  of  as  Pit  being 
the  matrix  indexed  by  (/+*’,  j+j’). 

Relative  pixel  intensity  between  the  two  images  can  be  used  as  a  measure  of 
match.  Alternatively,  correlation  is  well-suited  as  a  measure  of  match.  We  can 
mathematically  define  the  match  of  image  A  to  B,  Mab  ,  within  the  area  p,  as  the 
normalized  correlation  between  image  A  and  B,  or 


Mab(j,  j) 


2  X P0‘7)A(i  +  i\  j  +  /) B(i  + 1,  j  +  /) 


i  j  --m 


Sa(i,  j)  +  Ss(i,  j) 


,  [12] 


(47) 


where  m  can  range  from  0  to  2,  depending  upon  the  desired  area  of  p  (0  for  an  individual 
point,  2  for  a  5x5  matrix).  The  values  returned  for  match  are  between  -1  and  1.  Values 
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close  to  zero  indicate  low  correlation,  and  values  close  to  -1  or  1  indicate  high 
correlation. 

Each  level  of  the  pyramid  can  now  be  fused  by  testing  the  match  metric  between 
the  two  images  at  the  given  level.  If  the  match  metric  is  low  at  a  given  position,  then  the 
sample  from  the  source  image  with  the  highest  salience  is  copied  to  the  composite  image. 
If  the  match  metric  is  high,  the  pixels  from  the  source  images  are  averaged  and  copied  to 
the  composite  image.  This  combination  technique  can  be  implemented  as  a  weighted 
mean  in  which  the  weights  depend  upon  the  match  and  saliency  measures.  The  function 
FUSE,  as  described  above,  can  then  be  described  as  follows: 

C(i,  j)  =  um(i,  j)A(i,  j)  +  WB(i,  j)B(i,  j ) ,  (48) 

where  wa  and  wb  are  the  assigned  weights  to  the  source  images  A  and  B,  respectively; 
and  wa  +  wb  =  1 . 

There  are  two  common  ways  to  implement  the  weighting  functions  wa  and  wb  . 
One  way  is  to  select  some  threshold,  a,  for  the  match  metric.  We  can  set  wa  and  wb  to  1 
and  0  or  0  and  1,  respectively,  when  the  match  is  below  a,  or  we  can  set  wa  and  wb  to  .5 
and  .5  when  the  match  is  above  a.  This  would  be  mathematically  written  as 


if  a  <  threshold 
wa  =  .5  and  wb  =  .5, 
else  if  Sa  >  Sb 

wa  =  1  and  wb  —  0,  else  wa  =  0  and  wb  =  1. 


(49) 


This  technique  would  require  extensive  testing  to  ascertain  an  appropriate  a  for  each  type 
of  image  used  and  does  not  allow  for  a  gradual  change  from  selection  to  averaging. 
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Another  common  way  to  implement  the  weight  functions  is  to  assign  the  weight 
value  based  on  a  linear  transition  between  1  and  0  as  follows:  If  Mas  <  a ,  then  wmm  =  0 


and  Wmax  =  1 ;  otherwise, 


_  1 

Wmin  —  — 

2 


1  A-Mab^I 

2  1  -a 


and 


Wmax  —  1  Wmin  . 


(50) 

(51) 


WA  =  Wmax  and  W B  =  Wmin  if  Sa  >  Sb,  else  WA  =  Wmin  and  WB  =  Wmax  .  (52) 


Here,  the  larger  weight  is  assigned  to  the  source  image  with  the  higher  saliency  value. 
This  latter  approach  allows  a  gradual  change  from  selection  to  averaging  and  was  the 
method  used  in  this  study  because  it  offers  the  benefits  of  both  selection  and  averaging 
over  a  range  of  match  values  rather  than  simply  selecting  one  method  or  the  other.  This 
is  the  same  method  used  in  [12]. 
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CHAPTER  V 

FUSED  IMAGERY  AND  RESULTS 


A.  Fusion  Approach 

In  this  study,  six  multiresolution  image  fusion  functions  are  applied  using  the 
pyramid  and  fusion  techniques  described  in  Chapters  III  and  IV  to  demonstrate  the 
feasibility  of  image  fusion  in  remote  sensing  imagery.  The  first  three  functions  have 
been  described  in  the  current  literature.  The  other  three  functions  are  combinations  of 
pyramids  and  fusion  techniques  that  have  not  been  described  in  current  literature. 
Appendix  C  contains  computer  programs  used  to  implement  the  six  combinations  for 
multiresolution  image  fusion.  For  simplicity,  the  images  presented  in  this  chapter  are 
only  selected  subimages  of  the  images  presented  in  Chapter  II. 

Three  of  the  basic  pyramid  types  described  in  Chapter  III  for  image 
decomposition  are  employed  in  this  study.  The  pyramids  used  in  the  fusion  functions  are 
the  RoLP,  gradient,  and  morphological  pyramids.  (The  other  pyramids  in  Chapter  HI 
were  discussed  because  understanding  them  is  necessary  in  order  to  understand  the  three 
pyramids  that  were  used.)  These  three  pyramid  techniques  can  be  combined  with  the  two 
fusion  techniques  in  six  possible  ways.  The  three  combinations  described  in  the  current 
literature  are:  (1)  using  a  RoLP  pyramid  for  image  decomposition  with  contrast  fusion  for 
image  merging  [11,  16];  (2)  using  a  gradient  pyramid  for  image  decomposition  with 
hybrid  averaging  and  selection  fusion  for  image  merging  [12];  and  (3)  using  a 
morphological  pyramid  for  image  decomposition  with  contrast  fusion  for  image  merging 
[13].  The  three  not  discussed  in  current  literature  are:  (1)  using  a  RoLP  pyramid  for 
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image  decomposition  with  hybrid  averaging  and  selection  fusion  for  image  merging;  (2) 
using  a  gradient  pyramid  for  image  decomposition  with  contrast  fusion  for  image 
merging;  and  (3)  using  a  morphological  pyramid  for  image  decomposition  with  hybrid 
averaging  and  selection  fusion  for  image  merging. 

In  previous  studies,  composite  Landsat  and  SAR  images  have  been  created  by 
assigning  specific  bands,  or  ratios  of  bands,  to  specific  pixel  colors.  Images  that  have 
pixels  representing  both  SAR  and  Landsat  data  will  be  called  “hybrid  SAR  and  Landsat 
images”  or  simply  “hybrid  images.”  Hybrid  images  can  be  very  useful,  as  presented  in 
[5];  however,  they  limit  the  number  of  bands  that  can  be  viewed  to  three.  By  applying 
multiresolution  fusion  functions  to  source  images,  it  is  the  purpose  of  this  study  to  make 
it  possible  to  effectively  view  images  that  contain  information  from  more  than  three 
spectral  bands. 

The  remainder  of  this  thesis  presents  the  results  and  conclusions  of  applying  the 
six  fusion  functions  to  test  images.  Fusion  of  remotely  sensed  data  can  take  on  several 
forms.  In  the  following,  Section  B  presents  the  results  of  fusing  SAR  images  with 
Landsat  images  to  obtain  information  about  the  fusion  functions.  Section  C  discusses 
fusion  of  hybrid  images  to  obtain  useful  information  as  well  as  hybrid  image  fusion 
applications.  Finally,  in  Section  D,  the  factors  relevant  to  composite  image  quality  are 
considered. 

B.  Fusing  SAR  and  Landsat  Images  to 

Evaluate  Fusion  Functions 

In  this  section,  fusing  imagery  of  the  Mantua,  Utah,  area  is  reported.  This  set  of 
imagery  is  used  because  it  contains  mountainous  and  agricultural  areas  as  well  as  urban 
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areas.  A  list  of  the  figures  and  parameters  used  in  this  section  is  given  in  Table  III.  The 
first  source  image.  Figure  6  (a),  was  obtained  from  the  Landsat-4  satellite  on  25 
September  1994;  the  RGB  pixel  values  are  from  bands  4,  2,  and  1,  respectively.  The 
second  source  image,  Figure  6  (b),  was  obtained  from  the  SIR-C  SAR  instrument  flown 
on  STS-68  on  2  October  1994;  the  RGB  pixel  values  are  from  bands  L,  C,  and  C, 
respectively. 

The  results  presented  in  Figures  7,  8,  and  9  demonstrate  that  the  gradient  and 
RoLP  pyramids  with  the  hybrid  fusion  technique.  Figures  7  (b)  and  8  (b),  provide  the 
composite  images  with  the  most  well  integrated  fusion  of  source  image  features.  The 


TABLE  III 

LIST  OF  FIGURES  IN  CHAPTER  V,  SECTION  B 


Figure 

Area 

Description 

6(a) 

Mantua 

Landsat  imagery  of  Mantua,  bands  =  (4,  2, 1). 

6(b) 

Mantua 

SAR  imagery  of  Mantua,  bands  =  (L,  C,  C). 

7(a) 

Mantua 

Selection  fusion  of  6  (a)  and  6  (b)  using  a  four  level 
Gradient  pyramid  with  a  Gaussian  kernel  of  a  =  0.4. 

7(b) 

Mantua 

Hybrid  fusion  of  6  (a)  and  6  (b)  using  a  four  level 
Gradient  pyramid  with  a  Gaussian  kernel  of  a  =  0.4 
and  a  match  threshold  of  a  =  0. 1. 

8(a) 

Mantua 

Selection  fusion  of  6  (a)  and  6  (b)  using  a  four  level 
RoLP  pyramid  with  a  Gaussian  kernel  of  a  =  0.4. 

8(b) 

Mantua 

Hybrid  fusion  of  6  (a)  and  6  (b)  using  a  four  level 
RoLP  pyramid  with  a  Gaussian  kernel  of  a  =  0.4 
and  a  match  threshold  of  a  =  0. 1. 

9(a) 

Mantua 

Selection  fusion  of  6  (a)  and  6  (b)  using  a  four  level 
Morphological  pyramid  with  a  2x2  brick  as  a 
structuring  element. 

9(b) 

Mantua 

Hybrid  fusion  of  6  (a)  and  6  (b)  using  a  four  level 
Morphological  pyramid  with  a  match  threshold  of 
a  =  0. 1  and  a  2x2  brick  as  a  structuring  element . 

10  (a) 

Mantua 

Subsection  of  Figure  7  (b)  enlarged  for  detail. 

10(b) 

Mantua 

Subsection  of  Figure  8  (b)  enlarged  for  detail. 
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Fig.  6.  (a)  Landsat  and  (b)  SAR  source  images  of  the  Mantua,  Utah,  area. 
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(b) 

Fig.  7.  Fusion  of  the  Mantua,  Utah,  area  source  images.  Figure  6  (a)  and  6  (b),  using  a 
gradient  pyramid  and  (a)  selection  fusion,  (b)  hybrid  fusion. 
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visual  results  of  the  fusion  techniques  using  morphological  pyramids.  Figure  9,  show  that 
the  morphological  pyramids  do  not  work  well  for  remote  sensor  imagery.  The  low 
quality  of  the  composite  images  that  used  morphological  filters  is  due  to  the  small 
resolution  of  the  visual  elements;  for  example,  roads  are  only  one  or  two  pixels  wide.  In 
this  case,  using  a  2x2  or  a  3x3  structuring  element  for  the  morphological  filters  does  not 
allow  reconstruction  of  such  small  detail.  The  image  detail  cannot  be  reconstructed  with 
a  resolution  higher  than  that  of  the  structuring  element.  Hence,  the  composite  image  does 
not  contain  the  important,  small  details  from  the  source  images  because  the  details  are 
smaller  than  the  structuring  element.  A  smaller  structuring  element  is  not  feasible 
because  using  a  structuring  element  on  the  order  of  one  pixel  does  not  filter  the  image  at 
all.  Morphological  filters  work  well  for  high-resolution  images  where  image 
substructures  are  large  in  comparison  with  the  structuring  element.  However,  with  these 
images  the  entire  city  is  blurred. 

Results  obtained  using  the  selection  technique  for  image  fusion  are  unacceptable 
in  this  study  because  they  select  details  from  only  one  image  or  another.  In  the 
composite  images,  cities  and  mountains  end  up  being  represented  by  the  SAR 
information,  and  agricultural  areas  are  represented  by  the  Landsat  images.  This  is 
because  Landsat  images  have  the  highest  return  values  from  agricultural  areas,  and  SAR 
images  have  the  highest  return  values  in  cities  and  mountains.  Since  the  areas  of 
saturated  return  correspond  to  the  highest  pixel  values  and  the  contrast  selection 
technique  is  based  upon  ratios  of  pixel  intensity,  selection  fusion  leads  to  composite 
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images  that  have  the  SAR  information  in  city  areas  and  Landsat  information  in  the 
agricultural  areas. 

Figure  10  is  an  enlarged  subsection  of  the  two  fused  images  that  yielded  the  best 
result,  namely,  RoLP  and  gradient  pyramids  with  hybrid  fusion,  Figures  7  (b)  and  8  (b). 
Notice  that  the  detail  in  the  lower  left-hand  comer  of  the  RoLP  image  appears  blurry  as 
compared  with  the  same  area  in  the  gradient  composite  image. 

C.  Fusing  Hybrid  SAR/Landsat  Images 

In  the  previous  section,  SAR  and  Landsat  images  were  fused  to  obtain 
information  about  how  the  fusion  functions  work  for  remotely  sensed  imagery.  Fusion 
applications,  however,  are  more  likely  to  use  hybrid  images.  Hybrid  images  have  the 
advantage  of  using  multiple  sensors  to  obtain  single  source  images.  The  benefit  of 
having  the  different  bands  of  TM  sensors  and  SAR  sensors  is  that  each  band  shows  a 
particular  feature  of  the  surface.  When  these  features  are  understood,  source  images  can 
be  formed  to  show  specific  information  about  a  given  surface  area.  In  this  section,  the 
fusion  techniques  are  applied  to  a  sets  of  source  images  used  to  show  the  applicability  of 
fusion  of  hybrid  SAR  and  Landsat  images.  A  list  of  the  figures  used  in  this  section  is 
given  in  Table  IV. 

For  the  first  set  of  composite  images,  a  source  image  designed  to  show  urban, 
suburban,  and  agricultural  areas  is  combined  with  another  source  image  designed  to  show 
health  of  vegetation.  For  the  first  image,  we  want  to  use  a  band  that  reflects  vegetation 
and  a  band  that  reflects  anthropological  structures;  this  will  demonstrate  the  distinct 
difference  between  urban  and  agricultural  areas.  Landsat  band  2  is  in  the  visible 
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Fig.  10.  (a)  Gradient  pyramid  with  hybrid  fusion  enlarged  for  detail,  (b)  RoLP  pyramid 

with  hybrid  fusion  enlarged  for  detail. 


TABLE  IV 

LIST  OF  FIGURES  IN  CHAPTER  V,  SECTION  C 
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Figure 

Area 

Description 

11(a) 

Logan 

SAR  band  L  and  Landsat  bands  4  and  2. 

11(b) 

Logan 

Landsat  bands  3,  4,  and  2. 

12(a) 

Logan 

Hybrid  fusion  of  1 1  (a)  and  1 1  (b)  using  a  four  level 
Gradient  pyramid  with  a  Gaussian  kernel  of  a  =  0.4. 

12(b) 

Logan 

Hybrid  fusion  of  1 1  (a)  and  1 1  (b)  using  a  four  level 
RoLP  pyramid. 

spectrum  and  returns  a  peak  value  for  vegetation.  Landsat  band  4  is  in  the  near  infrared 
spectrum  and  shows  healthy  vegetation  and  land/water  interfaces.  Both  SAR  bands  C 
and  L  reflect  well  from  artificial  structures  and  would  work  well  for  this  image;  however, 
band  L  has  a  higher  return  from  artificial  structures  than  band  C,  and  band  C  reflects 
from  vegetation.  For  the  first  image,  we  use  Landsat  bands  2  and  4  and  SAR  band  L. 

The  second  source  image  is  designed  to  show  health  of  vegetation,  so  both  Landsat  bands 
2  and  4  are  again  used.  The  third  band  of  the  second  source  image  is  one  that  has  a  low 
return  in  vegetated  areas;  Landsat  band  3  is  in  the  visible-light  spectrum  that  corresponds 
to  chlorophyll  absorption. 

Figure  1 1  shows  the  two  hybrid  source  images.  The  first  image  emphasizes  land 
use  categories  and  land/water  boundaries,  and  the  second  shows  health  of  vegetation  and 
land/water  boundaries.  Figure  12  shows  the  result  of  fusing  the  source  images  using 
RoLP  and  gradient  pyramids  with  hybrid  fusion.  It  is  concluded  that  the  fusion  of  the 
source  images  is  successful  because  the  information  represented  by  Landsat  bands  2  and 
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49 


(b) 

Fig.  12.  Fusion  of  source  images  in  Figure  1 1  using  a  (a)  Gradient  pyramid  with  hybrid 
fusion,  (b)  RoLP  pyramid  with  hybrid  fusion. 
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4,  which  show  health  of  vegetation  and  water/land  boundaries,  remains  virtually 
unchanged,  whereas  the  composite  images  make  viewing  the  vegetation  in  the  city  easier 
without  changing  the  ease  of  viewing  the  land  use  information.  This  is  because  the  pixels 
represented  by  the  color  red  are  now  the  data  represented  by  the  fusion  between  the  SAR 
L  band  and  the  Landsat  band  3.  This  does  not  increase  the  amount  of  green  present  in  the 
image;  it  only  makes  it  easier  for  a  human  analyst  to  observe  because  of  the  way  we 
perceive  contrast.  The  reason  for  this  is  described  by  Weber’s  Law  [9].  By  decreasing 
the  amount  of  red  in  the  local  area,  it  decreases  the  amount  of  contrast  between  the  green 
and  red  pixels;  hence,  the  vegetation  (green)  is  more  easily  observed  by  a  human  analyst 
because  the  red  is  reduced.  If  we  were  to  simply  reduce  the  intensity  of  the  red  pixels, 
the  contrast  between  green  and  red  in  the  city  would  be  still  be  easily  observed;  however, 
in  areas  where  the  red  pixels  are  the  only  source,  the  contrast  would  also  be  reduced.  The 
use  of  fusion  allows  varying  the  amount  of  change  in  contrast  for  a  given  area  based  on 
the  correlation  between  the  two  source  images.  Figure  13  demonstrates  the  change  in 
contrast  by  viewing  enlarged  subsections  of  Figures  11  (a),  12  (a),  and  12  (b) 


(a)  (b)  (c) 

Fig.  13.  Enlarged  subsections  of  Figures  1 1  (a)  and  12  to  show  the  change  in  contrast. 
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respectively;  the  enlarged  subsections  in  Figure  13  are:  (a)  the  hybrid  image,  (b)  the  fused 
images  using  gradient  pyramids,  and  (c)  the  fused  images  using  RoLP  pyramids. 

To  most  observers  it  appears  that  Figure  13  (c)  has  the  most  green  and  Figure  13 
(a)  has  the  least  with  Figure  13  (b)  somewhere  in  between.  However,  all  three  images 
have  the  same  amount  of  green  in  them.  The  green  pixel  values  did  not  change;  only  the 
red  pixel  values  changed.  By  decreasing  the  intensity  of  the  red  pixel  values,  the  green  is 
more  easily  noticed;  hence,  people  see  more  green  in  the  images  (b)  and  (c)  of  Figure  13. 
The  same  effect  could  be  observed  over  small  areas  by  simply  decreasing  the  intensity  of 
the  red  pixels  in  the  image;  however,  fusion  allows  the  amount  of  change  in  the  pixel 
intensity  to  vary  dependent  upon  the  correlation  of  the  two  images.  Fusion  is  beneficial 
because  the  change  in  the  red  pixel  intensities  is  not  constant  over  the  whole  image.  The 
fusion  algorithm  decides  how  much  and  where  to  change  the  red  pixel  intensities. 

£>.  Composite  Image  Appearance  and  Quality 

The  composite  image  appearance  depends  upon  several  factors.  In  this  study, 
Landsat  and  SAR  images  were  converted  into  bitmaps  for  fusion.  This  preserved  the 
pixel  data.  In  order  for  the  REDUCE  and  EXPAND  functions  to  work  properly,  the  input 
image  size  needs  to  be  equal  to  a  power  of  two,  plus  one.  For  example,  the  two  sizes 
used  in  this  study  were  257x257  and  513x513.  This  works  because  257  =  28  +  1  and  513 
=  29  +  1.  If  a  source  image  is  passed  to  the  function  that  generates  a  pyramid,  the  image 
needs  to  have  x  and  y  sizes  that  are  a  power  of  two,  plus  one.  If  the  dimensions  of  the 
source  image  do  not  meet  this  criterion,  the  image  needs  to  be  resized.  When  resizing 
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occurs,  pixel  values  must  be  interpolated;  and  the  data  used  are  no  longer  exact. 
Therefore,  the  data  in  the  composite  image  are  not  exact. 

Another  consideration  in  the  composite  image  appearance  is  the  storage  format 
used.  If  an  image  is  stored  as  a  JPEG  file,  the  true  pixel  values  are  not  saved; 
quantization  is  necessary  for  the  compression  of  the  file.  Once  again,  because  the  data 
used  in  the  source  images  are  not  exact,  the  data  contained  in  the  composite  image  are  not 
exact.  It  is  necessary  for  the  user  to  decide  what  accuracy  of  pixel  values  is  necessary  for 
image  analysis  and  take  the  necessary  precautions  when  fusing  the  images. 

When  using  multiresolution  image  fusion,  the  number  of  levels  in  the  pyramids 
contributes  to  the  quality  of  the  composite  image.  For  example,  when  the  source  images 
are  decomposed  into  pyramids  six  levels  deep,  large  image  features  will  fuse  better  than 
if  pyramids  only  two  levels  deep  were  used.  The  depth  of  a  pyramid  is  an  important 
parameter.  If  the  depth  is  too  deep,  processing  time  is  wasted.  On  the  other  hand,  if  the 
pyramid  is  not  deep  enough,  the  larger  image  subfeatures  will  not  blend  well.  If  it  is 
known  that  the  remote  sensing  imagery  for  a  desired  study  has  surface/subsurface 
features  of  a  very  small  scale,  the  pyramid  depth  does  not  need  to  be  as  deep;  a  depth  of 
two  would  work  fine.  In  relation  to  the  overall  size  of  the  source  images  used  here,  the 
features  in  the  city  are  very  small,  while  the  features  in  the  mountains  are  quite  a  bit 
larger.  It  is  important  to  note  that  the  pyramid  depth  must  accommodate  the  fusion  of  the 
largest  subfeature  in  the  source  images;  this  is  what  will  determine  the  necessary  pyramid 
depth.  Figure  14  shows  images  created  by  fusing  two  images  using  gradient  pyramids  of 
different  depths.  Figure  15  shows  the  same  thing  using  RoLP  pyramids;  the  depths  used 
range  from  one  to  six. 
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Comparing  the  images  (in  Figure  14  for  gradient  pyramids  or  Figure  15  for  RoLP 
pyramids)  with  each  other,  it  can  be  seen  that,  after  a  pyramid  depth  of  two,  the  added 
pyramid  levels  do  not  visually  add  much  detail  to  the  composite  image.  The  most 
noticeable  change,  upon  visual  inspection,  is  in  the  color  of  the  mountains.  To  more 
clearly  see  the  change  that  results  from  increasing  the  pyramid  depth  from  level  to  level, 
Figure  16  shows  the  difference  in  gradient  images,  and  Figure  17  shows  the  difference  in 
RoLP  images.  (The  difference  of  each  pixel  value  has  been  multiplied  by  a  factor  of  4  to 
more  clearly  see  the  changes.) 

By  comparing  the  images  in  Figure  16  with  the  images  in  Figure  14  for  gradient 
pyramids  (or  Figure  17  with  Figure  15  for  RoLP  images),  it  is  noticed  that  the  difference 
from  level  to  level  is,  indeed,  in  the  mountainous  regions  of  the  image.  The  difference  in 
the  mountains  was  expected  because  the  image  features  in  the  mountains  are  much  larger 
than  the  image  features  in  the  city.  By  inspecting  the  resultant  images,  it  can  be  observed 
that  the  detail  in  the  composite  images  for  the  city  areas  did  not  change  noticeably  for  any 
of  the  pyramid  depths  used  after  a  depth  of  two.  For  any  set  of  images,  the  optimum 
pyramid  depth  depends  upon  the  size  of  the  details  considered  for  analysis.  The  larger 
the  details  in  the  source  images,  the  deeper  the  pyramid  depth  needs  to  be  for  satisfactory 
fusion.  If  the  pyramids  were  skipped  all  together  (a  pyramid  of  depth  =  1),  the  composite 
image  has  a  higher  level  of  detail  missing  as  opposed  to  using  a  depth  of  two. 

This  study  also  used  a  subjective  measure  of  image  quality  as  described  by  [24].  The 
quality  of  each  image  has  been  assigned  a  value  between  0  and  9.  An  image  with  a 
quality  measure  (QM)  of  1  would  be  able  to  locate  only  large  terrain  features  such  as  rail 
yards,  airstrips,  or  possibly  very  large  aircraft  at  an  airstrip.  An  image  with  a  QM  of  8, 
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(e) 

Fig.  16.  Four  times  the  difference  between  gradient  composite  images  of  successively 
increasing  pyramid  depths,  (a)  Fig.  14  (b)  -  Fig.  14  (a),  (b)  Fig.  14  (b)  -  Fig.  14  (c).  (c) 
Fig.  14  (d)  -  Fig.  14  (c).  (d)  Fig.  14  (e)  -  Fig.  14  (d).  (e)  Fig.  14  (f)  -  Fig.  14  (e). 
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(e) 

Fig.  17.  Four  times  the  difference  between  RoLP  composite  images  of  successively 
increasing  pyramid  depths,  (a)  Fig.  15  (b)  -  Fig.  15  (a),  (b)  Fig.  15  (b)  -  Fig.  15  (c).  (c) 
Fig.  15  (d)  -  Fig.  15  (c).  (d)  Fig.  15  (e)  -  Fig.  15  (d).  (e)  Fig.  15  (f)  -  Fig.  15  (e). 
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on  the  other  hand,  would  be  able  to  detect  types  of  aircraft  at  an  airstrip  or  the  number  of 
boxcars  at  a  rail  yard.  The  image  quality  was  obtained  from  the  power  spectral  density  of 
the  image  luminance.  The  image  luminance  is  acquired  by  transforming  from  the  RGB 
coordinate  system  to  the  XYZ  coordinate  system,  where  Y  is  the  image  luminance. 
(Appendix  B  contains  the  transform  equations  and  the  image  QMs  for  the  composite 
images  in  this  study.)  This  quality  measure  indicated  that  image  quality  decreased  as  the 
pyramid  depth  increased.  It  also  indicated  that  the  composite  morphological  images  had 
higher  quality  than  the  RoLP  and  gradient  composite  images,  which  is  not  the  case,  as 
can  be  visually  verified.  It  is,  therefore,  concluded  that  the  quality  measure  as  described 
in  [24]  does  not  work  well  for  fused  imagery  of  remotely  sensed  data. 
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CHAPTER  VI 
CONCLUSIONS 

A.  Findings 

The  findings  from  this  study  are  verified  by  the  fused  images  generated  from  SAR 
and  TM  data  of  the  same  terrestrial  scenes.  It  has  been  shown  that  there  is  a  more 
effective  way  to  view  composite  SAR  and  Landsat  images  than  by  simply  viewing  three 
bands  at  one  time.  The  cost  of  the  more  effective  composite  image  is  computation  time. 
Table  V  is  a  summary  of  processor  time  for  the  various  fusion  techniques  that  were 
studied. 

The  fastest  fusion  techniques  use  morphological  filters;  however,  as  the 
composite  images  indicate,  morphological  filters  do  not  work  well  for  these  types  of 
remotely  sensed  data  because  of  the  small-scale  of  image  subfeatures  in  comparison  with 
the  structuring  element. 


TABLE  V 

COMPUTATION  TIMES  FOR  FUSION  TECHNIQUES 

Technique 

Size  (Pixels) 

Time  (Seconds) 

Morphological-Selection 

513x513 

116.2 

RoLP-Selection 

513x513 

130.1 

Morphological-Hybrid 

513x513 

260.2 

RoLP-Hybrid 

513x513 

631.2 

Gradient-Selection 

513x513 

1142.6 

Gradient-Hybrid 

513x513 

1823.4 

Computational  time  on  a  Pentium  266  MHz  processor  with  64  MB  of  RAM. 
For  an  image  size  of  257x257,  all  computational  times  should  be  quartered. 
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Fusion  techniques  using  RoLP  pyramids  are  the  second  fastest.  The  composite 
images  formed  from  RoLP  techniques  are  beneficial  in  that  they  present  the  data  from 
both  source  images  in  the  composite  scene  better  than  do  the  simple  hybrid  images  used 
in  previous  studies.  Compared  with  techniques  that  use  gradient  pyramids,  however,  the 
composite  images  from  RoLP  pyramids  tend  to  appear  blurry.  This  appearance  of  being 
blurry,  however,  is  an  artifact  of  the  reduced  contrast  in  the  composite  image.  The 
quality  remains  approximately  constant  between  the  gradient  and  RoLP  composite 
images  (given  equal  pyramid  depths). 

Gradient  pyramids  yielded  composite  images  in  approximately  equal  quality  to,  or 
slightly  better  than,  the  RoLP  pyramids  (gradient  composites  do  not  appear  as  blurry); 
however,  the  computation  time  for  gradient  pyramids  is  about  four  times  that  of  the  RoLP 
pyramids.  The  gradient  pyramid  is  actually  a  combination  of  four  pyramids.  When  exact 
analysis  of  imagery  is  necessary,  the  gradient  pyramid  may  be  the  best  choice;  however, 
for  most  applications,  the  RoLP  pyramid  with  hybrid  fusion  will  work  just  as  well. 

It  was  found  in  this  study  that  the  selection  approach  to  image  fusion  is 
undesirable.  Because  the  selection  technique  chooses  only  one  band  or  another  to 
represent,  it  does  not  fuse  them.  Although  the  selection  technique  is  almost  twice  as  fast 
as  the  hybrid  selection  and  averaging  technique,  the  hybrid  approach  yields  a  composite 
image  with  better  detail.  A  result  of  this  study,  therefore,  is  the  suggestion  that,  for 
general  fusion  applications,  a  RoLP  pyramid  with  hybrid  selection  and  averaging  fusion 
should  be  used.  If  small  details  are  a  concern  and  time  is  not  a  constraint,  the  gradient 
pyramid  with  hybrid  fusion  may  yield  a  slightly  better  result. 
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This  study  has  demonstrated  that  the  optimum  pyramid  depth  depends  upon  the 
largest  important  subfeature  in  the  image.  If  the  features  of  concern  are  small  in  detail, 
like  the  cities  in  the  case  of  satellite  data,  a  pyramid  depth  of  two  will  produce  the  same 
result  as  a  pyramid  depth  of  six.  On  the  other  hand,  if  the  features  in  the  image  are  large, 
like  mountains,  a  much  deeper  pyramid  is  necessary. 

A  simple  rule  to  follow  would  be  to  use  a  pyramid  depth  of  two  if  the  image 
features  desired  for  fusion  are  on  the  order  of  one  to  two  pixels.  For  most  applications,  a 
pyramid  depth  of  three  would  work  fine.  If  the  image  features  are  very  large,  on  the 
order  of  hundreds  of  pixels,  a  pyramid  depth  of  six  would  be  appropriate.  Rarely  would  a 
pyramid  depth  of  greater  than  five  or  six  be  needed.  Based  upon  the  results  from  figures 
16  and  17,  a  pyramid  depth  of  at  least  two  should  be  used  in  all  fusion  applications 
because  the  amount  of  detail  added  when  changing  from  a  pyramid  of  depth  one  to  a 
pyramid  of  depth  two  is  sufficiently  large. 

B.  Recommendations 

Future  research  related  to  the  topic  of  this  study  is  needed  in  several  areas.  First, 
an  optimum  alpha  value  (the  threshold  that  is  used  to  regulate  averaging  and  selection  of 
composite  pixel  values)  for  the  match  and  saliency  measure  could  be  established  through 
research.  New  approaches  need  to  be  developed  to  speed  up  the  gradient  filters,  thus 
greatly  benefiting  the  software  already  developed.  Also,  blind  source  separation  [25]  can 
be  explored  as  a  means  to  image  fusion.  Another  area  of  research  is  to  find  more  suitable 
structuring  elements  for  morphological  filters.  Structuring  elements  that  would  make 
possible  the  retention  of  small-scale  image  features  in  the  composite  image  are  needed. 
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Another  area  of  a  study  could  be  the  comparison  between  using  the  wavelet 
transform  or  using  pyramids  for  multiresolution  representation  of  images.  Another  type 
of  pyramid/fusion  combination  that  could  be  tried  is  a  ratio  between  levels  of  a 
morphologically  filtered  pyramid  with  both  fusion  types.  The  final  suggestion  I  have  for 
further  research  is  to  develop  another  way  to  objectively  measure  image  quality  without 
knowing  the  image  source.  Such  a  metric  would  be  a  valuable  tool  for  the  field  of  image 
processing  and  analysis  in  general.  Further  research  in  these  areas  could  greatly  benefit 
the  field  of  image  processing  as  well  as  the  analysis  of  remote  sensing  imagery. 
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Knowing  that,  theoretically,  the  RoLP  pyramid  should  have  an  exact 
reconstruction  and  that  morphological  and  gradient  pyramids  should  have  approximate 
reconstructions,  the  following  tests  were  performed:  First,  a  small  test  image  of  17x17 
pixels  was  created;  then  the  pyramid  routines  were  applied;  and,  at  each  stage  of  the 
image  decomposition  and  reconstruction,  the  values  were  verified  by  hand  calculations. 
The  second  test  was  to  take  three  larger  test  images,  shown  in  Figure  18,  create  a 
pyramid,  and  then  recover  the  image  by  using  the  pyramid  reconstruction  techniques. 
Finally,  the  reconstructed  image  was  subtracted  from  the  original  image;  and  the  average 
and  standard  deviation  of  pixel  differences  was  calculated.  These  values  are  given  for  the 
test  images  in  Figure  1 8  in  Tables  VI,  VII,  and  VIII.  The  second  test  was  run  with 
pyramids  of  depth  =  6  and  pyramids  of  depth  =  4.  For  an  exact  reconstruction,  the 
average  pixel  difference  would  be  zero;  and  the  standard  deviation  would  also  be  zero. 
The  testing  showed  that  the  RoLP  pyramid  was  exact  in  most  cases.  The  one  exception 
was  for  a  six  level  pyramid  for  image  18  (b)  in  the  red  band;  in  which  case  it  was  still 
very  near  zero.  See  Table  VII.  The  morpological  pyramids  had  the  largest  errors,  and 
the  errors  were  independent  of  pyramid  depth.  The  errors  in  morphological  pyramids  are 
due  to  the  fact  that  detail  smaller  than  the  structuring  element  cannot  be  exactly 
reconstructed.  Gradient  pyramids  had  a  moderate  margin  of  error  for  pyramid 
reconstruction;  these  errors  are  due  to  the  approximation  shown  in  equations  (19)  and 
(20).  In  general,  the  deeper  the  pyramid  the  larger  the  error  in  reconstructing  the  image. 
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(a)  (b)  (c) 


Fig.  18.  Test  images  used  to  verify  pyramid  code  correctness. 


TABLE  VI 

RESULTS  FROM  TESTING  PYRAMID 
RECONSTRUCTION  CODE  ON  FIGURE  19  (a) 


Pyramid  Type  Pyramid  Depth  Pixels  Color 


TABLE  VII 

RESULTS  FROM  TESTING  PYRAMID 
RECONSTRUCTION  CODE  ON  FIGURE  19  (b) 
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Pyramid  Type 

Pyramid  Depth 

Pixels  Color 

Average  Error  Standard  Deviation 

Morph 

6 

Red 

6.048186 

16.756359 

Morph 

6 

Green 

2.045051 

4.711925 

Morph 

6 

Blue 

0.810453 

1.817184 

RoLP 

6 

Red 

0.000418 

0.051939 

RoLP 

6 

Green 

0.000000 

0.000000 

RoLP 

6 

Blue 

0.000000 

0.000000 

Gradient 

6 

Red 

4.250710 

4.461939 

Gradient 

6 

Green 

5.615083 

5.245300 

Gradient 

6 

Blue 

4.639384 

1.951420 

Morph 

4 

Red 

6.048186 

16.756359 

Morph 

4 

Green 

2.045051 

4.711925 

Morph 

4 

Blue 

0.810453 

1.817184 

RoLP 

4 

Red 

0.000418 

0.051939 

RoLP 

4 

Green 

0.000000 

0.000000 

RoLP 

4 

Blue 

0.000000 

0.000000 

Gradient 

4 

Red 

2.729921 

2.905167 

Gradient 

4 

Green 

3.289907 

2.882259 

Gradient 

4 

Blue 

2.969148 

1.107655 

TABLE  VIII 

RESULTS  FROM  TESTING  PYRAMID 

RECONSTRUCTION  CODE  ON  FIGURE  19  (c) 

Pyramid  Type 

Pyramid  Depth 

Pixels  Color 

Average  Error  Standard  Deviation 

Morph 

6 

Red 

1.887326 

5.636188 

Morph 

6 

Green 

1.660585 

4.812446 

Morph 

6 

Blue 

1.723327 

4.875687 

RoLP 

6 

Red 

0.000000 

0.000000 

RoLP 

6 

Green 

0.000000 

0.000000 

RoLP 

6 

Blue 

0.000000 

0.000000 

Gradient 

6 

Red 

5.413866 

6.271135 

Gradient 

6 

Green 

4.584484 

3.704715 

Gradient 

6 

Blue 

4.457820 

3.544401 

Morph 

4 

Red 

1.887326 

5.636188 

Morph 

4 

Green 

1.660585 

4.812446 

Morph 

4 

Blue 

1.723327 

4.875687 

RoLP 

4 

Red 

0.000000 

0.000000 

RoLP 

4 

Green 

0.000000 

0.000000 

RoLP 

4 

Blue 

0.000000 

0.000000 

Gradient 

4 

Red 

3.306019 

3.456579 

Gradient 

4 

Green 

3.044988 

2.178596 

Gradient 

4 

Blue 

3.016866 

2.076221 

Appendix  B:  Image  Quality  Measures 
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The  image  quality  measures  used  in  this  study  were  calculated  from  the  power 

spectrum  of  the  image  as  described  in  [24].  The  power  spectrum  was  taken  on  the  image 

luminance,  which  was  calculated  by  transforming  from  the  RGB  coordinates  to  the  XYZ 

coordinates  through  the  following  transform: 

X  =  0.490*R  +  0.310*G  +  0.200*B; 

Y  =  0.177*R  +  0.813*G  +  0.01 1*B; 

Z  =  0.000*R  +  0.010*G  +  0.990*B; 

where  Y  is  the  image  luminance. 

The  quality  measure  was  found  as  follows: 

QM  =  1.93*log  (IQM)  +  8.77, 

where 

1  180°  0.5 

IQM  =  £  JjS(dl)W(p)A2 (Tp)P(p,d)  [24]. 

M  0=180°  P=0.01 

Here,  P(p,0)  is  the  Discrete  Fourier  Transform  or  the  Fast  Fourier  Transform.  5(0,)  is 
the  scale  of  the  image  at  an  angle  off?,  from  the  sensor.  W  (p)  is  the  classic  Wiener  filter 
as  described  in  [24];  and  A(Tp )  is  a  model  of  the  human  visual  system  (HVS),  also 
described  in  [24], 

The  quality  measures  (QMs)  were  taken  by  dividing  images  into  quadrants  of  size 
128x128.  An  image  of  original  size  513x513  would  be  divided  into  quadrants  and  then 
subdivided  again  into  subquadrants.  Quadrants  start  with  number  one  in  the  upper  right- 
hand  comer  of  the  image  and  then  proceed  counterclockwise  to  the  lower  right-hand 
comer,  which  is  quadrant  four.  Table  IX  and  X  list  the  QMs  for  the  composite  images 
created  in  the  study.  The  QMs  from  this  study  indicate  that  a  quantitative  measure  of 
image  quality  that  models  the  HVS  exactly  still  needs  to  be  developed.  This  quality 
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measure  did,  however,  provide  a  means  to  measure  small  differences  in  image  quality 
and  trends  in  changes  based  upon  varying  inputs  such  as  pyramid  depth. 


TABLE  IX 

QUALITY  MEASURES  FOR  COMPOSITE  IMAGES  IN  FIGURES  7  THRU  12 


Figure 

Quadrant 

S_Q1 

S_Q2 

S_Q3 

S_Q4 

7(a) 

i 

5.743940 

6.299867 

5.955395 

5.798092 

7(a) 

2 

5.774915 

6.010545 

6.144174 

6.240824 

7(a) 

3 

6.031102 

6.109628 

5.829595 

6.149844 

7(a) 

4 

5.537493 

4.949940 

5.178320 

6.203316 

7(b) 

1 

5.731352 

5.888458 

5.774751 

5.797175 

7(b) 

2 

5.455630 

5.956551 

5.839167 

6.143242 

7(b) 

3 

5.871299 

5.776086 

6.154200 

6.210767 

7(b) 

4 

5.686476 

5.879713 

5.157750 

6.044588 

8(a) 

1 

6.216465 

6.068850 

6.105093 

5.860565 

8(a) 

2 

2.987027 

5.877355 

5.704545 

6.401773 

8(a) 

3 

6.054478 

6.100063 

6.158357 

6.205588 

8(a) 

4 

3.824745 

4.298860 

5.471602 

6.395981 

8(b) 

1 

5.683046 

6.041419 

5.824888 

5.818359 

8(b) 

2 

5.043132 

5.987381 

5.779756 

6.235353 

8(b) 

3 

5.879769 

5.643103 

6.231497 

6.232010 

8(b) 

4 

5.597918 

5.909996 

5.156698 

6.147477 

9(a) 

1 

5.542208 

6.417702 

5.934047 

5.910607 

9(a) 

2 

6.183458 

5.925222 

4.363320 

6.441932 

9(a) 

3 

6.170932 

6.095793 

6.064970 

6.194372 

9(a) 

4 

4.960927 

5.389463 

5.162599 

6.263077 

9(b) 

1 

5.757378 

6.148040 

5.833893 

5.915841 

9(b) 

2 

6.171905 

6.070212 

5.708892 

6.401994 

9(b) 

3 

5.996606 

5.932110 

6.246206 

6.177732 

9(b) 

4 

4.377702 

5.238811 

5.299767 

6.259186 

11(a) 

1 

4.182329 

6.433823 

6.187418 

6.323740 

11(a) 

2 

5.828491 

5.019885 

5.546538 

5.774281 

11(a) 

3 

5.602613 

6.041763 

5.954587 

4.099891 

11(a) 

4 

6.281577 

6.463826 

5.614375 

5.843056 

11(b) 

1 

4.182329 

6.433823 

6.187418 

6.323740 

11(b) 

2 

5.828491 

5.019885 

5.546538 

5.774281 

11(b) 

3 

5.602613 

6.041763 

5.954587 

4.099891 

11(b) 

4 

6.281577 

6.463826 

5.614375 

5.843056 

12(a) 

1 

4.971820 

6.396582 

6.191291 

6.352974 

12(a) 

2 

5.787445 

4.995888 

5.494674 

5.742789 

12(a) 

3 

5.592110 

5.999828 

5.919660 

4.059719 

12(a) 

4 

6.216224 

6.469916 

5.645947 

5.817126 

12(b) 

1 

4.170645 

6.430707 

6.190845 

6.319584 

12(b) 

2 

5.826273 

5.018929 

5.545628 

5.773333 

12(b) 

3 

5.604957 

6.040998 

5.953582 

4.104306 

12(b) 

4 

6.282585 

6.463218 

5.615763 

5.842026 

Here  S_Q1  stands  for  subquadrant  I,  S_Q2  stands  for  subquadrant  2,  S_Q3  stands  for 
_ subquadrant  3,  and  S_Q4  stands  for  subquadrant  4. _ 
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TABLE  X 

QUALITY  MEASURES  FOR  COMPOSITE  IMAGES  IN  FIGURES  14  AND  15 


Figure 

Quadrant 

S_Q1 

S_Q2 

S_Q3 

S_Q4 

14(a) 

i 

5.907000 

6.105254 

5.954466 

5.892723 

14(a) 

2 

6.134328 

6.087080 

6.248683 

6.403800 

14(a) 

3 

6.060922 

5.777319 

6.310244 

6.218068 

14(a) 

4 

5.394370 

5.815682 

5.238525 

6.211088 

14(b) 

1 

5.909584 

6.089230 

5.919499 

5.878862 

14(b) 

2 

6.041452 

6.063083 

6.129162 

6.358926 

14(b) 

3 

6.009208 

5.737093 

6.263228 

6.211189 

14(b) 

4 

5.493500 

5.796942 

5.213909 

6.161926 

14(c) 

1 

5.860167 

6.027700 

5.870211 

5.852189 

14(c) 

2 

5.782127 

6.023773 

5.988911 

6.284763 

14  (c) 

3 

5.948137 

5.751993 

6.220279 

6.218034 

14  (c) 

4 

5.596609 

5.785432 

5.199620 

6.111913 

14(d) 

1 

5.731352 

5.888458 

5.774751 

5.797175 

14(d) 

2 

5.455630 

5.956551 

5.839167 

6.143242 

14  (d) 

3 

5.871299 

5.776086 

6.154200 

6.210767 

14(d) 

4 

5.686476 

5.879713 

5.157750 

6.044588 

14(e) 

1 

5.449706 

5.570606 

5.653515 

5.714482 

14  (e) 

2 

5.736781 

5.890406 

5.472869 

5.945735 

14  (e) 

3 

5.797911 

5.752635 

6.020828 

6.168142 

14  (e) 

4 

5.714855 

6.015885 

5.034786 

5.996046 

14(f) 

1 

5.155452 

5.912589 

5.555240 

5.590639 

14(f) 

2 

6.083034 

5.791876 

4.703432 

5.677327 

14(f) 

3 

5.674115 

5.736203 

5.816588 

6.045265 

14(f) 

4 

5.674343 

6.074334 

4.679347 

5.855797 

15(a) 

1 

5.907000 

6.105254 

5.954466 

5.892723 

15(a) 

2 

6.134328 

6.087080 

6.248683 

6.403800 

15(a) 

3 

6.060922 

5.777319 

6.310244 

6.218068 

15(a) 

4 

5.394370 

5.815682 

5.238525 

6.211088 

15(b) 

1 

5.894845 

6.126315 

5.929002 

5.887159 

15(b) 

2 

5.856588 

6.067307 

6.129728 

6.374213 

15(b) 

3 

6.004285 

5.757083 

6.264523 

6.222906 

15(b) 

4 

5.463844 

5.787212 

5.208783 

6.179084 

15(c) 

1 

5.840663 

6.041954 

5.889577 

5.861843 

15(c) 

2 

5.188877 

6.028841 

5.963924 

6.309322 

15(c) 

3 

5.939251 

5.734792 

6.252913 

6.221472 

15(c) 

4 

5.520334 

5.778367 

5.183959 

6.152245 

15(d) 

1 

5.683046 

6.041419 

5.824888 

5.818359 

15(d) 

2 

5.043132 

5.987381 

5.779756 

6.235353 

15  (d) 

3 

5.879769 

5.643103 

6.231497 

6.232010 

15(d) 

4 

5.597918 

5.909996 

5.156698 

6.147477 

15(e) 

1 

5.313802 

5.692202 

5.792126 

5.804777 

15(e) 

2 

5.652391 

5.971659 

5.431761 

6.216579 

15(e) 

3 

5.864790 

5.557136 

6.161402 

6.245972 

15(e) 

4 

5.673966 

6.072573 

5.161863 

6.130646 

15(f) 

1 

4.772618 

5.091612 

5.772559 

5.813745 

15(f) 

2 

6.126851 

5.980894 

4.989271 

6.207485 

15(0 

3 

5.839331 

5.599754 

6.058500 

6.202442 

15(f) 

4 

5.702899 

6.123158 

5.102581 

6.125611 

Here  S_Q1  stands  for  subquadrant  1,  S_Q2  stands  for  subquadrant  2,  S_Q3  stands  for 
subquadrant  3,  and  S_Q4  stands  for  subquadrant  4. 

Appendix  C:  Source  Code 


/********************************★*★  fuse_n. cpp  ********************************/ 
/*  This  file  contains  the  code  for  fusing  the  images.  It  contains  both  */ 

/*  the  fusion  functions  for  selection  and  hybrid  selection  and  averaging.  */ 

/*  It  also  contains  the  six  fusion  functions  used  to  fuse  RoLP,  gradient,  */ 

/*  and  morphological  pyramids  with  the  two  fusion  techniques.  */ 

/*******************************************************************************/ 
#ifndef  FUSE 
#def ine  FUSE 

#include  " thesis. h* 
tinclude  "files.cpp" 

#include  "convolve_n. cppw 
tinclude  "gradient_n. cpp" 

# include  *morph_n.cpp" 

# include  *rolp_n. cpp" 


/*******************************************************************************/ 
/*  Image  Fusion  Techniques  */ 

/*************♦*****************************************************************/ 
/*  Fuses  two  images  selecting  between  the  two  source  pixels  the  pixel  with 
the  highest  intensity;  also  called  contrast  fusion.  */ 
image  *  selection (image  *il ,  image  *i2) 

{ 

image  *newpic; 

int  xnew,  ynew,  x,  y,  index; 
if (il  ==  NULL  ||  i2  ==  NULL)  return (NULL) ; 
if (il->xsize  ==  i2->xsize  &&  il“>ysize  ==  i2->ysize)  { 
xnew  =  il->xsize; 
ynew  =  i2->ysize; 

newpic  =  (image  *)malloc (sizeof (image) ) ; 

newpic->pixel  =  (Pixel  *)malloc (xnew*ynew* sizeof (Pixel) ) ; 
for(x  =  0;  x  <  xnew;  x++)  { 

for(y  =  0;  y  <  ynew;  y++)  { 

index  =  y*xnew  +  x; 

if(  ( il->pixel [index] [0]  -  1)  >  ( i2->pixel [index] [0 ]  -  1))  { 

newpic->pixel [index] [0 ]  =  il->pixel [index] [0 ] ; 

} 

else  { 

newpic ->pixel [index]  [0]  =  i2->pixel [index]  [0 ]  ; 

} 

if (  (il->pixel [index] [1]  -  1)  >  (i2->pixel [index] [1]  -1))  { 

newpic ->pixel [index] [1]  =  il->pixel [index] [1] ; 

} 

else  { 

newpic ->pixel [index] [1]  =  i2->pixel [index] [1] ; 

} 

if (  (il->pixel [index] [2]  -  1)  >  (i2->pixel [index] [2]  -  1))  { 

newpic“>pixel [index] [2]  =  il->pixel [index] [2 ] ; 

} 

else  { 

newpic->pixel [index] [2]  =  i2->pixel [index] [2] ; 

} 

} 

} 

newpic->xsize  =  xnew; 
newpic->ysize  =  ynew; 
return (newpic) ; 

> 

else  return (NULL); 


/*  This  function  is  used  for  the  match  and  saliency  measures  for  hybrid 
selection  and  averaging  fusion.  This  function  returns  the  weighted 
average  of  the  product  of  a  3 -by- 3  area  surrounding  the  pixel  in  the 
(if  j)  position  for  a  given  pixel  color.  */ 
double  s amp le_p (image  *  picl,  image  *  pic2,  int  i,  int  j,  int  color) 

{ 

//The  images  must  be  the  same  x  and  y  dimentions 


int  m,  n,  index,  tx,  ty; 
double  tsum  =  0.0; 
for(m  =  -2;  m  <  3;  m++)  { 

for{n  =  -2;  n  <  3;  n++)  { 

tx  =  i  +  m; 

if (tx  <0  ||  tx  >=  picl->xsize  | |  tx  >=  pic2->xsize)  continue; 
ty  =  j  +  n; 

if ( ty  <0  ||  ty  >=  picl->ysize  ||  ty  >=  pic2->ysize)  continue; 
index  =  ty*picl->xsize  +  tx; 

tsum  +=  w_g(m,n)  *  picl->pixel [index] [color]  *  pic2->pixel [index] [color] ; 

} 

} 

return  (tsum)  ; 

} 

/*  The  salience  function  returns  a  corresponding  map  the  same  dimentions  of  the 
image  which  is  used  to  determine  which  pixel  to  use  for  selection,  or  which 
one  has  a  higher  weight  for  averaging.  */ 
image  *  salience (image  *g) 

{ 

image  *  newp i c ; 
int  x,  y,  index; 

newpic  =  (image  * )malloc (sizeof ( image) ) ; 

newpic->pixel  =  (Pixel  *)malloc (g->xsize*g->ysize*sizeof (Pixel) ) ; 
newpic->xsize  =  g->xsize; 
newpic->ysize  =  g->ysize; 
for (y  =  0;  y  <  g->ysize;  y++)  { 

for (x  =  0;  x  <  g->xsize;  x++)  { 

index  =  y*g->xsize  +  x; 

newpic->pixel [index] [0]  =  sample_p(g,  g,  x,  y,  0) ; 

newpic->pixel [index] [1]  =  sample_p(g,  g,  x,  y,  1); 

newpic->pixel [index] [2 ]  =  sample_p(g,  g,  x,  y,  2); 

> 

} 

return ( newp ic ) ? 

} 

/*  This  function  reuturns  a  measure  of  match  to  by  which  to  compare  the  two 
images.  If  the  two  images  match,  averaging  is  used,  if  they  don’t  match 
within  the  weighted  average  is  used  where  the  one  with  the  higher 
saliency  has  the  higher  weight.  */ 
image  *  match (image  *a,  image  *b,  image  *sa,  image  *sb) 

{ 

image  *newpic; 
int  x,  y,  index; 

newpic  =  (image  *)malloc (sizeof (image) ) ; 

newpic->pixel  =  (Pixel  * Jmalloc (a->xsize*a->ysize*sizeof (Pixel )) ; 
newpic->xsize  =  a->xsize; 
newpic->ysize  =  a->ysize; 
for(y  =  0;  y  <  a->ysize;  y++)  { 

for(x  =  0;  x  <  a->xsize;  x++)  { 

index  =  y*a->xsize  +  x; 

newpic->pixel [ index]  [0 ]  =  abs (2*sample_p  (a,  b,  x,  y,  0)/ 

(sa->pixel [index] [0]  +  sb->pixel [index] [ 0 ] ) ) ; 
newpic->pixel [index]  [1]  =  abs (2*sample_p (a,  b,  x,  y,  1)/ 

(sa->pixel [index]  [1]  +  sb->pixel [index]  [  1 ] ) ) ; 
newpic- >pixel [index] [2]  =  abs (2*sample_p  (a,  b,  x,  y,  2)/ 

(sa->pixel [index]  [2]  +  sb->pixel [index]  [2] ) ) ; 

> 

} 

return (newpic) ; 

} 

/*  Fuses  two  images  by  using  both  selection  and  averaging.  If  the  two  pixels 
correspond  to  the  same  type  of  subfeature  averageing  is  used,  otherwise 


selection  is  used.  Alpha  is  a  measure  of  how  close  the  two  pixels  must 
match  to  use  selection  or  averaging.  */ 
image  *  hybrid__selection_aver age  (image  *a,  image  *b,  double  alpha) 

{ 

image  *sa,  *sb,  *m,  *newpic; 
double  wmin,  wmax,  wa,  wb? 
int  x,  y,  index,  c; 

if (a  ==  NULL  | |  b  ==  NULL)  return (NULL) ; 
sa  =  salience (a); 
sb  =  salience (b) ; 
m  =  match(a,  b,  sa,  sb)  ; 

newpic  =  (image  *)malloc (sizeof (image) ) ; 

newpic->pixel  =  (Pixel  *)malloc (a->xsize*a->ysize*sizeof (Pixel) ) ; 
newpic->xsize  =  a->xsize; 
newpic->ysize  =  a->ysize; 
for (y  =  0;  y  <  a->ysize;  y++)  { 

for(x  =  0;  x  <  a->xsize;  x++)  { 

index  =  y*a->xsize  +  x; 
for(c  =  0;  c  <  3;  C++)  { 

wmin  =  .5  -  .5*(1  -  m->pixel [index] [c] )/ (1  -  alpha); 
wmax  =  1  -  wmin; 

if (sa->pixel [index] [c]  >  sb->pixel [index] [c] )  { 

wa  =  wmax;  wb  =  wmin; 

}  else  { 

wa  =  wmin;  wb  =  wmax; 

} 

newpic ->pixel [index] [c]  =  wa*a->pixel [index] [c]  +  wb*b->pixel [index] [c] ; 

} 

} 

} 


return (newpic) ; 

} 

/*******************************************************************************/ 

/*  Fuse  Gradient  Pyramid  with  hybrid  fusion  */ 

/★★★★★★★a***********************************************************************/ 

/*  This  function  fuses  two  gradient  pyramids  using  hybrid  fusion  */ 

image  *  fuse_gradient_hybrid__pyramid  (Pyramid_Gradient  *A,  Pyramid_Gradient  *B,  double 
alpha) 

{ 

Pyramid„Gradient  *C; 
image  *comp; 


C  =  (Pyramid_Gradient  * )malloc (sizeof (Pyramid_Gradient )) ; 
C->d4  =  (Pyramid  *)malloc (sizeof (Pyramid) ) ; 

C->d3  =  (Pyramid  * Jmalloc (sizeof (Pyramid) ) ; 

C->d2  =  (Pyramid  * )malloc (sizeof (Pyramid) ) ; 

C->dl  =  (Pyramid  *)malloc (sizeof (Pyramid) ) ; 

C->g  =  (Pyramid  *)malloc (sizeof (Pyramid) ) ; 


print f ( " \nFusing  Gradient  4 . . . ■ ) ; 

C->d4->I5  =  NULL; 

C->d4->I4  =  hybrid_selection_average (A->d4->I4, 
C->d4->l3  =  hybrid_selection_average (A->d4->I3 , 
C->d4->l2  =  hybrid_selection_average(A->d4->I2, 
C~>d4->I1  =  hybrid_selection_average (A->d4->Il, 
C->d4->l0  =  hybrid_selection_average (A->d4->I0 , 
print f ( " . . .Done\n" ) ; 


B->d4->I4 , 
B->d4->I3 , 
B->d4->I2, 
B->d4->Il , 
B->d4->I0 , 


alpha)  ; 
alpha)  ; 
alpha) ; 
alpha) ; 
alpha)  ; 


print f ( " \nFusing  Gradient  3  .  .  .  " )  ; 

C->d3->I5  =  NULL; 

C->d3->l4  =  hybrid_selection_average (A->d3->I4, 
C->d3->I3  =  hybrid_selection_average (A->d3->I3 , 
C->d3->I2  =  hybrid_selection_average (A->d3->I2 , 
C->d3->Il  =  hybrid_selection_average (A->d3->Il, 
C->d3->I0  =  hybridise lection_average (A->d3->I0, 
printf ( ■ . . .Done\n" ) ; 


B->d3->I4 , 
B->d3->I3 , 
B->d3->I2 , 
B->d3->Il , 
B->d3->I0, 


alpha)  ; 
alpha) ; 
alpha)  ; 
alpha)  ; 
alpha) ; 


printf ( " \nFusing  Gradient  2  .  .  .  " ) ; 

C->d2->I5  =  NULL; 

C->d2->l4  =  hybrid_selection_average  (A->d2-*>I4, 
C->d2->l3  =  hybrid_selection_average  (A->d2*->I3 , 
C->d2->I2  =  hybrid_selection_average (A->d2~>I2, 
C->d2->Il  =  hybrid_selection_average (A->d2->Il, 
C->d2->I0  =  hybrid_selection_average (A->d2->lO , 
printf ( " . . .DoneXn" ) ; 


B->d2->I4, 
B->d2->I3/ 
B->d2->I2 , 
B->d2->Il , 
B->d2->!0, 


alpha) ; 
alpha)  ; 
alpha)  ; 
alpha)  ; 
alpha)  ; 


printf ( * \nFusing  Gradient  1 . .  .  " )  ; 

C->dl->I5  =  NULL; 

C->dl->l4  =  hybrid__selection_average  (A->dl->l4,  B->dl->l4, 
C->dl->I3  =  hybrid_selection_average (A->dl->I3 ,  B->dl->I3, 
C->dl->I2  =  hybrid_selection__average  (A->dl->I2 ,  B->dl->l2 , 

C->dl->Il  =  hybrid_selection_average (A->dl->Il,  B->dl->Il, 
C”>dl->10  =  hybrid_selection_average (A->dl->I0/  B->dl->I0, 
printf ( " . . .Done\n" ) ; 


alpha) ; 
alpha) ; 
alpha) ; 
alpha) ; 
alpha) ; 


printf ( " \nFusing  Gradient  0  .  . .  " ) ? 

C->g->I5  =  hybrid_selection_average (A->g->I5/  B->g->I5/  alpha); 
C->g->l4  =  NULL; 

C->g->I3  =  NULL; 

C->g->I2  =  NULL; 

C->g->Il  =  NULL; 

C->g->lO  =  NULL; 
printf ( * . . . Done\n" ) ; 


printf { " \nConstructing  Composite  Pyramid. . . " ) ; 
comp  =  const ruct_gradi ent (C) ; 
free_Pyramid_Gradient (C) ; 
printf ( " . . .Done\n" ) ; 


return (comp) ; 

} 

/*  This  function  fuses  two  images  by  constructing  gradient  pyramids  and 
then  fusing  the  pyramids  using  hybrid  fusion.  */ 
image  *  fuse_gradient_hybrid( image  *il,  image  *i2,  double  alpha) 

{ 

image  *c; 

Pyramid_Gradient  *A,  *B; 

A  =  (Pyramid_Gradient  * )malloc (sizeof (Pyramid_Gradient ) ) ; 

B  =  (Pyramid_Gradient  *)malloc (sizeof (Pyramid_Gradient) ) ; 

printf ( * \nCreating  Pyramid  A. . . " ) ; 

A  =  gen _j>yramid_gradient ( il ) ; 

printf ( " . . .DoneNn" ) ; 

printf ( • XnCreating  Pyramid  B . . . " ) ; 

B  =  gen  pyramid  gradient (i2) ; 
printf ( * . . . DoneXn* ) ; 

c  =  fuse_gradient_hybrid_pyramid(A,  B,  alpha); 

f ree_Pyramid_Gradient  (A)  ; 
free_Pyramid_Gradient (B) ; 


return ( c ) ; 

} 

/★★★★★★★★a**********************************************************************/ 

/*  Fuse  RoLp  Pyramid  with  selection  fusion  */ 

/a******************************************************************************/ 

/*  This  function  fuses  two  RoLP  pyramids  using  selection  fusion  */ 
image  *  fuse_RoLP_selection_pyramid (Pyramid  *A,  Pyramid  *B) 

{ 

Pyramid  *C; 
image  *comp; 


81 


C  =  (Pyramid  *)malloc (sizeof (Pyramid) ) ; 

C->I5  =  selection (A->I5,  B->I5) ; 

C->I4  =  selection (A->I4,  B->I4); 

C->I3  =  selection (A->I3 ,  B->I3); 

C->I2  =  selection (A->I2,  B->I2); 

C->I1  =  selection (A->I1,  B->I1)  ; 

C->I0  =  selection (A->I0#  B->I0); 

comp  =  cons t rue tJRoLP (C) ; 
free_pyramid(C) ; 

return (comp) ; 

} 

/*  This  function  fuses  two  images  by  constructing  RoLP  pyramids  and  then 
fusing  the  pyramids  using  selection  fusion.  */ 
image  *  fuse_RoLP_se lection (image  *il,  image  *12 ) 

{ 

image  *c; 

Pyramid  *A,  *B; 

A  =  (Pyramid  *)malloc (sizeof (Pyramid) ) ; 

B  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

A  =  gen_pyramid_RoLP (il)  ; 

B  =  gen_pyramid_RoLP ( 12 ) ; 

c  =  fuse__RoLP__selection_pyramid  (A,  B)  ; 

free_pyramid(A) ? 
free_pyramid(B)  ; 

return (c) ; 

} 

/**********************************************************************■*********/ 
/*  Fuse  Morph  Pyramid  with  selection  fusion  */ 

y*******************************************************************************/ 

/*  This  function  fuses  two  morphological  pyramids  using  selection  fusion  */ 
image  *  fuse„morph_selection_pyramid (Pyramid  *A,  Pyramid  *B) 

{ 

Pyramid  *C ; 
image  *comp; 

C  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

C->I5  =  selection (A->I5/  B->I5) ; 

C->I4  =  selection (A->I4 ,  B->I4) ; 

C-’>I3  =  selection  (A->I3/  B->I3); 

C->I2  =  selection (A->I2 ,  B->I2); 

C->I1  =  selection (A->I1/  B->I1) ; 

C->I0  =  selection (A->I0 7  B->I0); 

comp  =  construct„morph(C); 
free_pyramid(C) ; 

return (comp) ; 

} 

/*  This  function  fuses  two  images  by  constructing  morphological  pyramids  and 
then  fusing  the  pyramids  using  selection  fusion.  */ 
image  *  fuse_morph_select ion ( image  *il,  image  *i2) 

{ 

image  *c; 

Pyramid  *A,  *B; 

A  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

B  =  (Pyramid  *)malloc (sizeof (Pyramid)); 


init_sK() ; 
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A  =  gen_pyramid_morph(il) ; 

B  =  gen_pyramid_morph  ( i2 )  ; 

c  =  fuse_morph_selection_pyramid(A,  B) ; 

free  ^pyramid  (A)  ; 
free_pyramid(B)  ; 


return ( c ) ; 

} 

/*******************************************************************************/ 
/*  Fuse  Morph  Pyramid  with  hybrid  fusion  */ 

/********************************************************************** *********/ 
/*  This  function  fuses  two  morphological  pyramids  using  hybrid  fusion  */ 
image  *  fuse_morph_hybrid_pyramid( Pyramid  *A,  Pyramid  *B,  double  alpha) 

{ 

Pyramid  *C; 
image  *comp; 


C  =  (Pyramid  *) malloc (sizeof (Pyramid)); 


C->I5  =  hybrid_selection_average (A->I5, 
C->I4  =  hybrid_selection_average (A->I4 , 
C->I3  =  hybrid_selection_average (A->I3 , 
C->I2  =  hybrid_selection_average  (A->I2  , 
C->I1  =  hybrid_selection_average (A->I1/ 
C->I0  =  hybrid_selection_average (A->I0 , 


B->I5 ,  alpha); 
B->I4,  alpha); 
B~>13 ,  alpha); 
B->I2,  alpha); 
B->I1 ,  alpha); 
B->I0 ,  alpha); 


comp  =  construct_morph(C); 
free_pyramid(C)  ; 

return (comp) ; 

} 


/*  This  function  fuses  two  images  by  constructing  morphological  pyramids  and 
then  fusing  the  pyramids  using  hybrid  fusion.  */ 
image  *  fuse_morph_hybr id ( image  *il,  image  *i2,  double  alpha) 

{ 

image  *c; 

Pyramid  *A,  *B; 

A  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

B  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 


init_sK ( ) ; 

A  =  gen_pyramid_morph(il) ; 

B  =  gen_pyramid_morph ( i2 ) ; 

c  =  fuse_morph_hybrid_pyramid(A,  B,  alpha); 


free  ^pyramid  (A)  ; 
free_pyramid(B)  ; 


return ( c ) ; 

} 

/*******************************************************************************/ 
/*  Fuse  RoLp  Pyramid  with  hybrid  fusion  */ 

/*******************************************************************************/ 
/*  This  function  fuses  two  RoLP  pyramids  using  hybrid  fusion  */ 
image  *  fuse_RoLP„hybrid_pyramid( Pyramid  *A,  Pyramid  *B#  double  alpha) 

{ 

Pyramid  *C; 
image  *comp; 

C  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

C->I5  =  hybrid_selection_average (A~>15/  B->I5,  alpha); 

C->I4  =  hybrid__selection_average  (A->I4,  B->I4,  alpha); 


C->I3 

0>I2 

C->I1 

C->IO 


hybrid_selection_average (A->I3 ,  B->I3 , 
hybrid_selection_average ( A->I2 ,  B->I2 , 
hybrid_selection„average ( A->I1 ,  B~>I1 , 
hybrid_selection_average {A- >10 ,  B->I0 , 


alpha)  ? 
alpha)  ; 
alpha)  ; 
alpha)  ; 


comp  =  cons t rue t_RoLP (C) ; 
free  _pyramid(C) ; 


return (comp) ; 

} 

/*  This  function  fuses  two  images  by  constructing  RoLP  pyramids  and 
then  fusing  the  pyramids  using  hybrid  fusion.  */ 
image  *  fuse_RoLP_hybrid (image  *il,  image  *i2,  double  alpha) 

{ 

image  *c; 

Pyramid  *A,  *B; 

A  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

B  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

A  =  gen_pyramid_RoLP (il) ; 

B  =  gen_pyramid__RoLP  (i2)  ; 

c  =  f  us  e_RoLP_hybr  id_pyramid ( A ,  B,  a lpha ) ? 

free_pyramid(A) ; 
free_pyramid(B) ; 

return ( c ) ; 

> 

/**************************************************★****************************/ 
/*  Fuse  Gradient  Pyramid  with  hybrid  fusion  */ 

/a******************************************************************************/ 
/*  This  function  fuses  two  gradient  pyramids  using  selection  fusion  */ 
image  *  fuse_gradient__selection__pyramid  (Pyramid_Gradient  *A,  Pyramid„Gradient  *B) 
{ 

Pyramid_Gradient  *C ; 
image  *comp; 

C  =  (Pyramid__Gradient  *)malloc  (sizeof  (Pyramid_Gradient)  )  ; 

C->d4  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

C->d3  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

C->d2  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

C->dl  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

C->g  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 


C->d4->I5 

C->d4->I4 

C->d4->I3 

C->d4->I2 

C->d4->Il 

C->d4->I0 


NULL; 

selection ( A->d4->l4 , 
selection ( A->d4->l3 , 
selection (A->d4->I2 , 
selection (A->d4->Il , 
selection (A->d4->I0 , 


B->d4->l4 ) 
B->d4->l3 ) 
B->d4->l2 ) 
B->d4->Il) 
B->d4->!0 ) 


C->d3->I5 

C~>d3->I4 

C->d3->l3 

C->d3->l2 

C->d3->Il 

C->d3->I0 


NULL; 

selection (A->d3->I4, 
selection (A->d3->I3 , 
selection (A->d3->I2 , 
selection (A->d3->Il, 
selection (A->d3->I0, 


B->d3->I4) ; 
B->d3->I3 )  ; 
B->d3->I2)  ; 
B->d3->Il)  ; 
B->d3->!0 )  ; 


C->d2->l5 

C->d2->l4 

C->d2->I3 

C->d2->I2 

C->d2->Il 

C->d2->I0 


NULL; 

selection ( A->d2->I4 , 
selection (A->d2->I3 , 
selection (A->d2->I2 , 
selection (A->d2->Il, 
selection ( A->d2->!0 , 


B->d2->I4 ) 
B->d2->I3 ) 
B->d2->l2) 
B->d2->Il ) 
B->d2->!0) 


C->dl->I5  =  NULL; 
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C->dl->I4  =  selection (A->dl->l4,  B->dl->I4) ; 

C->dl->I3  =  selection (A->dl->I3 ,  B->dl->I3); 

C->dl->l2  =  selection (A->dl->l2 ,  B->dl->l2); 

C->dl->Il  =  selection (A->dl->Il,  B->dl->Il) ; 

C->dl->lO  =  selection (A->dl->I0 ,  B->dl->I0); 

C->g->I5  =  selection (A->g->I5#  B->g->I5) ; 

C->g->l4  =  NULL; 

C->g->I3  =  NULL; 

C->g->I2  =  NULL; 

C->g->H  =  NULL; 

C->g->lO  =  NULL; 

comp  =  construct_gradient(C); 
free_Pyramid„Gradient (C) ; 

return (comp) ; 

} 

/*  This  function  fuses  two  images  by  constructing  gradient  pyramids  and 
then  fusing  the  pyramids  using  selection  fusion.  */ 
image  *  fuse_gradient_selection { image  *il,  image  *i2) 

{ 

image  *c; 

Pyramid__Gradient  *A,  *B; 

A  =  (Pyramid_Gradient  * )malloc (sizeof (Pyramid_Gradient) ) ; 

B  =  (Pyramid_Gradient  *)malloc (sizeof (Pyramid_Gradient) ) ; 

print f ( " \nCreating  Pyramid  A. . . " ) ; 

A  =  gen_pyramid_gradient (il) ; 

printf ( ’ . . .Done\nw ) ; 

print f ( " \nCreating  Pyramid  B . . . " ) ; 

B  =  gen_pyramid_gradient (i2 ) ; 
printf ( * . . .Done\n" ) ; 

c  =  fuse_gradient_selection_jpyramid(A/  B)  ? 

free_Pyramid_Gradient (A) ; 
free_Pyramid_Gradient (B) ; 

return ( c ) ? 

} 


#endif 
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/************************************  image_p. cpp  ******************************/ 
/*  This  file  contains  the  code  used  for  image  processing  routines  such  as  */ 
/*  image  addition,  subtraction,  multiplication,  division,  etc.  */ 

/♦it*****************************************************************************/ 

#ifndef  IMAGE 
#define  IMAGE 

# include  "thesis.h" 

#include  'files. cpp* 

/*  Type  Definitions  */ 
typedef  struct  Pyramid  { 

image  *10,  *11,  *12,  *13,  *14,  *15; 

}  Pyramid; 

/*  Image  Processing  */ 

image  *  reduce (image  *g) ; 

image  *  expand (image  *g) ; 

image  *  subtract ( image  *il,  image  *i2); 

image  *  add (image  *il,  image  *i2) ; 

image  *  max_comp ( image  *il,  image  *i2); 

image  *  divide (image  *il,  image  *i2); 

image  *  multiply ( image  *il,  image  *i2) ; 

void  free„pyramid( Pyramid  *p)  ; 

/★★★a***************************************************************************/ 
/*  IMAGE  PROCESSING  STUFF  */ 

/A******************************************************************************/ 

/*  The  Gaussian  kernel,  w.  */ 
double  w(int  m,  int  n) 

{ 

double  wm,  wn,  answer; 

if (m  ==  0)  wm  =  A„REDUCE ; 

else  if  (abs(m)  ==  1)  wm  =  0.25; 

else  wm  =  0.25  -  A_REDUCE/2; 

if (n  ==  0)  wn  =  A_REDUCE ; 

else  if  (abs(n)  ==  1)  wn  =  0.25; 

else  wn  =  0.25  -  A_REDUCE/2; 

answer  =  wm*wn; 
return (answer) ; 

} 

/*  Returns  a  single  sample  located  at  (i,  j)  for  the  REDUCED  image  of  pic 
color  is  the  indicator  for  red,  green,  or  blue  pixels.  */ 
int  reducepixel (image  *  pic,  int  i,  int  j,  int  color) 

{ 

int  m,  n,  index,  tx,  ty; 
double  tsum  =  0.0; 
for(m  =  -2;  m  <  3;  m++)  { 

for(n  =  -2;  n  <  3;  n++)  { 

tx  =  2*i  +  m; 
if(tx  <  0)  tx  =  -tx; 

if(tx  >=  pic->xsize)  tx  =  2*pic->xsize  -  tx  -  1; 

ty  =  2* j  +  n; 

if (ty  <  0)  ty  =  -ty; 

if(ty  >=  pic->ysize)  ty  =  2*pic->ysize  -  ty  -  1; 

index  =  ty*pic->xsize  +  tx; 

tsum  +=  w(m,n)  *  pic->pixel [index] [color] ; 

} 

} 

return ( (int) tsum) ; 

} 

/*  Filters  and  subsamples  the  image  g  to  get  the  next  level  in  the 
Gaussian  pyramid.  */ 


image  *  reduce (image  *g) 

{ 

image  *  newp i c ; 

int  xnew,  ynew,  x,  y,  index; 

newpic  =  (image  *)malloc (sizeof (image) ) ; 

xnew  =  (g->xsize  -  1 ) / 2  +  1; 

ynew  =  (g->ysize  -  l)/2  +  1; 

newpic->pixel  =  (Pixel  * )malloc (xnew*ynew* sizeof (Pixel) ) ; 
for(x  =  0;  x  <  xnew;  x++)  { 

for (y  =  0;  y  <  ynew;  y++)  { 

index  =  y*xnew  +  x; 

newpic->pixel [index] [0]  =  reducepixel (g,  x,  y,  0); 
newpic->pixel [index] [1]  =  reducepixel (g,  x,  y,  1); 
newpic ->pixel [index] [2]  =  reducepixel (g,  x,  y#  2); 

} 

} 

newpic->xsize  =  xnew; 
newpic ->ysize  =  ynew; 
return (newpic) ; 


/*  Returns  a  single  sample  located  at  (i,  j)  for  the  EXPANDED  image  of  pic 
color  is  the  indicator  for  red,  green,  or  blue  pixels.  */ 
double  expandpixel (image  *  pic,  int  i,  int  j,  int  color) 

{ 

int  m,  n,  index,  tx,  ty; 
double  tsum  =  0.0; 
for(m  =  -2;  m  <  3;  m++)  { 

for(n  =  -2;  n  <  3;  n++)  { 

if (  ( (i-m)  %2  ==  0)  &&  ((j-n)%2  ==  0)  )  { 

tx  =  (i-m)  /  2; 
if(tx  <  0)  tx  =  -tx; 

if ( tx  >=  pic->xsize)  tx  =  2*pic->xsize  -  tx  -  1; 

ty  =  (j  -  n)/2; 
if(ty  <  0)  ty  =  -ty; 

if(ty  >=  pic->ysize)  ty  =  2*pic->ysize  -  ty  -  1; 

index  =  ty*pic->xsize  +  tx; 

tsum  +=  w(m,n)  *  pic->pixel [index] [color] ; 

> 

} 

} 

return (4*tsum) ; 


/*  Upsamples  and  interpolates  the  image  g  to  get  a  lower  resolution  image 
of  the  previoius  level  in  the  pyramid*/ 
image  *  expand (image  *g) 

{ 

image  * newpic; 

int  xnew,  ynew,  x,  y,  index; 

newpic  =  (image  *)malloc (sizeof (image) ) ; 

xnew  =  (g->xsize  -  1)*2  +  1; 

ynew  =  (g->ysize  -  1)*2  +  1; 

newpic->pixel  =  (Pixel  *) malloc (xnew*ynew*sizeof (Pixel )) ; 
for(x  =  0;  x  <  xnew;  x++)  { 

for(y  =  0;  y  <  ynew;  y++)  { 

index  =  y*xnew  +  x; 

newpic->pixel [index] [0]  =  expandpixel (g,  x,  y,  0); 
newpic->pixel [index] [1]  =  expandpixel (g,  x,  y,  1); 
newpic->pixel [index] [2]  =  expandpixel (g,  x,  y,  2); 

} 

} 

newpic->xsize  =  xnew; 
newpic->ysize  =  ynew; 
return (newpic) ; 
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/*  This  function  returns  the  absolut  value  of  the  difference  between 
image  il  and  image  i2 .  */ 

image  *  subtract (image  *il ,  image  *i2) 

{ 

image  *newpic; 

int  xnew,  ynew,  x,  y,  index; 

if (il->xsize  ==  i2->xsize  &&  il->ysize  ==  i2~>ysize)  { 
xnew  =  il->xsize; 
ynew  =  i2->ysize; 

newpic  =  (image  * )malloc (sizeof (image) ) ; 

newpic->pixel  =  (Pixel  *) malloc (xnew*ynew* sizeof (Pixel) ) ; 
for(x  =  0;  x  <  xnew;  x++)  { 

for(y  =  0;  y  <  ynew;  y++)  { 

index  =  y*xnew  +  x; 

newpic->pixel [index] [0]  =  abs (il->pixel [index] [0]  -  i2->pixel [index] [0 ]) ; 
newpic->pixel [index] [1]  =  abs (il->pixel [index] [1]  -  i2->pixel [index] [1] ) ; 
newpic ->pixel [index] [2]  =  abs (il->pixel [index] [2]  -  i2->pixel [index] [2] ) ; 

} 

} 

newpic->xsize  =  xnew; 
newpic ->ysize  =  ynew; 
return (newpic) ; 

} 

else  return (NULL); 

> 

/*  This  function  returns  the  difference  between  image  il  and  image  i2 .  */ 

image  *  subtract_sign ( image  *il,  image  *i2) 

{ 

image  *  newp i c ; 

int  xnew,  ynew,  x,  y,  index; 

if(il->xsize  ==  i2->xsize  &&  il->ysize  ==  i2->ysize)  { 
xnew  =  il->xsize; 
ynew  =  i2->ysize; 

newpic  =  (image  *) malloc (sizeof ( image) ) ; 

newpic->pixel  =  (Pixel  *) malloc (xnew*ynew* sizeof (Pixel) ) ; 
for(x  =  0;  x  <  xnew;  x++)  { 

for (y  =  0;  y  <  ynew;  y++)  { 

index  =  y*xnew  +x; 

newp ic->pixel [index] [0]  =  il->pixel [index] [0 ]  -  i2->pixel [index] [0] ; 
newpic- >pixel [index] [1]  =  il->pixel [index] [1]  -  i2->pixel[index][l]; 
newpic- >pixel [index] [2]  =  il->pixel [index] [2 ]  -  i2->pixel [ index] [2] ; 

> 

} 

newpic->xsize  =  xnew; 
newpic->ysize  =  ynew; 
return (newpic) ; 

} 

else  return (NULL) ? 

} 

/*  This  function  returns  the  sum  of  image  il  and  image  i2.  */ 

image  *  add_sign (image  *il,  image  *i2) 

{ 

image  *newpic; 

int  xnew,  ynew,  x,  y,  index; 

if(il->xsize  ==  i2->xsize  &&  il->ysize  ==  i2->ysize)  { 
xnew  =  il->xsize; 
ynew  =  i2->ysize; 

newpic  =  (image  *) malloc (sizeof ( image) ) ; 

newpic->pixel  =  (Pixel  *) malloc (xnew*ynew*sizeof (Pixel) ) ; 
for(y  =  0;  y  <  ynew;  y++)  { 

for(x  =  0;  x  <  xnew;  x++)  { 

index  =  y*xnew  +  x; 

newpic->pixel [index] [0 ]  =  il->pixel [index] [0 ]  +  i2->pixel [index] [0] ; 
newpic->pixel [index] [1]  =  il->pixel [index] [ 1]  +  i2->pixel [index] [1] ; 
newpic->pixel [index] [2]  =  il->pixel [index] [2]  +  i2->pixel [index] [2] ; 


} 


} 

newpic ->xsize  =  xnew; 
newpic->ysize  =  ynew; 
return (newpic) ; 

} 

else  return (NULL) ; 

} 

/*  This  function  frees  the  memory  used  by  a  pyramid.  */ 
void  f ree_pyramid ( Pyramid  *p) 

{ 

if (p  ! =  NULL)  { 
free_image (p->10 ) ; 
free_image (p->Il) ; 
free_image (p->l2) ; 
free_image (p->l3 ) ; 
free_image (p->14) ; 
free_image (p->15) ; 
free (p) ; 
p  =  NULL; 

} 

return ; 

} 

/*  This  function  returns  the  ratio  of  image  il  and  image  i2 .  */ 

image  *  divide (image  *il,  image  *i2) 

{ 

image  *newpic ; 

int  xnew,  ynew,  x,  y,  index; 

if (il->xsize  ==  i2->xsize  &&  il->ysize  ==  i2->ysize)  { 
xnew  =  il->xsize; 
ynew  =  i2->ysize; 

newpic  =  (image  * )malloc (sizeof (image) ) ; 

newpic->pixel  =  (Pixel  *)malloc (xnew*ynew*sizeof (Pixel) ) ; 
for(x  =  0;  x  <  xnew;  x++)  { 

for(y  =  0;  y  <  ynew;  y++)  { 

index  =  y*xnew  +  x; 
if (i2->pixel [index] [0]  !=  0)  { 

newpic->pixel [index] [0]  =  (il->pixel [index] [0]  /  i2->pixel [index] [0 ] ) 

} 

else  { 

newpic->pixel [index] [0]  =  1; 

} 

if (i2->pixel [index] [1]  !=  0)  { 

newpic->pixel [index] [1]  =  ( il->pixel [index] [1]  /  i2->pixel [index] [1] ) 

} 

else  { 

newpic ->pixel [index] [1]  =  1; 

} 

if (i2->pixel [index] [2]  !=  0)  { 

newpic ->pixel [index] [2]  =  (il->pixel [index] [2 ]  /  i2->pixel [index] [2] ) 

} 

else  { 

newpic- >pixel [index] [2]  =  1; 

} 

} 

} 

newpic->xsize  =  xnew; 
newpic->ysize  =  ynew; 
return (newpic) ; 

} 

else  return (NULL) ; 

} 

/*  This  function  returns  the  product  of  image  il  and  image  i2 .  */ 

image  *  multiply (image  *il,  image  *i2) 

{ 


image  * newpic; 

int  xnew,  ynew,  x,  y,  index; 

if(il->xsize  ==  i2->xsize  &&  il->ysize  ==  i2->ysize)  { 
xnew  =  il->xsize; 
ynew  =  i2->ysize; 

newpic  =  {image  *) malloc {sizeof (image) ) ; 

newpic->pixel  =  (Pixel  *) malloc (xnew*ynew*sizeof (Pixel) ) ; 
for (x  =  0;  x  <  xnew;  x++)  { 

for(y  =  0;  y  <  ynew;  y++)  { 

index  =  y*xnew  +  x; 

newpic->pixel  [index]  [0]  =  (il->pixel  [index]  [0]  *  i2-*>pixel  [index]  [0 ]  ) 
newpic->pixel [index] [1]  =  (il->pixel [index] [1]  *  i2->pixel [index] [1] ) 
newpic->pixel [index] [2]  =  (il->pixel [index] [2]  *  i2->pixel [index] [2] ) 

} 

} 

newpic->xsize  =  xnew; 
newpic->ysize  =  ynew; 
return (newpic) ; 

} 

else  return (NULL) ; 

} 

/*  This  function  creates  an  image  structure  from  the  bands  of  other  images, 
il,  i2,  and  i3  are  the  images  used,  and  cl,  c2,  and  c3  are  the  bands  that 
will  form  the  new  red,  green,  and  blue  bands.  */ 
image  *  spec_bands ( image  *il,  int  cl,  image  *i2,  int  c2,  image  *i3,  int  c3) 

{ 

image  * newpic; 

int  xnew,  ynew,  x,  y,  index; 

if ( il->xsize  ==  i2->xsize  &&  il->ysize  ==  i2->ysize  && 
i2->xsize  ==  i3->xsize  &&  i2->ysize  ==  i3->ysize)  { 
xnew  =  il->xsize; 
ynew  =  i2~>ysize; 

newpic  =  (image  *) malloc (sizeof (image) ) ; 

newpic->pixel  =  (Pixel  *) malloc (xnew*ynew*sizeof (Pixel) ) ; 
for(x  =  0;  x  <  xnew;  x++)  { 

for(y  =  0;  y  <  ynew;  y++)  { 

index  =  y*xnew  +  x; 

newpic->pixel [index] [0]  =  il->pixel [index] [cl] ; 
newpic->pixel [index] [1]  =  i2->pixel [index] [c2] ; 
newpic->pixel [index] [2]  =  i3->pixel [index] [c3] ; 

} 

} 

newpic- >xsize  =  xnew; 
newpic->ysize  =  ynew; 
return (newpic) ; 

> 

else  return(NULL) ; 


/*  This  function  creates  a  ppm  file  from  the  bands  of  other  ppm  files. 
The  new  image  contains  the  red  band  of  nl,  the  green  band  of  n2,  and 
the  blue  band  of  n3.  */ 

void  mergebands ( char  *  out,  char  *  nl,  char  *  n2,  char  *  n3) 

{ 

int  xsize,  ysize,  x,  y; 

FILE  *  infilel,  *infile2,  *infile3,  *outfile; 
int  red,  green,  blue,  tl,  t2; 
initfile (&infilel,  nl,  &xsize,  &ysize); 
initfile (&infile2,  n2,  &xsize,  &ysize); 
initfile (&infile3 ,  n3 ,  &xsize,  &ysize) ; 
outfile  =  fopen(out,  "w* ) ; 

fprintf (outf ile,  "P3\n#  Created  by  Ted  Meek\n"); 
fprintf (outfile,  "%d  %d  %d\n",  xsize,  ysize,  255); 
for (y  =  0;  y  <  ysize;  y++)  { 

for(x  =  0;  x  <  xsize;  x++)  { 

fscanf (infilel,  ■%d%d%d",  &red,  &tl,  &t2); 
f scanf { inf ile2 ,  "%d%d%d",  &tl,  &green,  &t2); 
fscanf (inf ile3,  "%d%d%d",  &tl,  &t2,  &blue) ; 
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fprintf (outfile,  "%d  %d  %d  ",  red,  green,  blue); 

} 

} 

fclose (inf ilel) ; 
f close ( inf ile2 ) ? 
fclose (inf ile3) ; 
fclose (outfile) ; 

} 

/*  This  function  returns  the  product  of  image  il  and  the  constand  c.  */ 
image  *  multiply_const (image  *il ,  double  c) 

{ 

image  *newpic; 

int  xnew,  ynew,  x,  y,  index; 
xnew  =  il->xsize; 
ynew  =  il->ysize; 

newpic  =  (image  *) malloc (sizeof (image) ) ; 

newpic->pixel  =  (Pixel  * )malloc (xnew*ynew* sizeof (Pixel) ) ; 
for(x  =  0;  x  <  xnew;  x++)  { 

for(y  =  0;  y  <  ynew;  y++)  { 

index  =  y*xnew  +  x; 

newpic->pixel [index] [0]  =  (il->pixel [index] [0]  *  c); 
newpic->pixel [index] [1]  =  (il->pixel [index] [1]  *  c); 
newpic->pixel [index] [2]  =  (il->pixel [index] [2]  *  c) ; 
for (int  i  =  0;  i  <3;  i++)  { 

if ( newpic- >pixel [index] [i]  >  255)  newpic->pixel [index] [i]  =  255; 

} 

} 

} 

newpic->xsize  =  xnew; 
newpic->ysize  =  ynew; 
return (newpic) ; 


#endif 
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/**********************************  gauss ian_n . cpp  *****************************/ 

/*  This  file  contains  the  code  for  generating  a  Gaussian  pyramid  for  a  given  */ 

/*  image.  The  pyramid  depth  is  six  levels.  */ 

/★★★a***************************************************************************/ 

#ifndef  GAUSSIAN 
#def ine  GAUSSIAN 

# include  " thesis. h" 

#include  "files. cpp" 

# include  " image_p . cpp " 

/***********************★**********************************************************/ 
/*  Gaussian  */ 

/******************************************★*************************************★*/ 
/*  This  functin  generates  a  six  level  Gaussian  pyramid  for  the  given  image  */ 

Pyramid  *  gen_pyramid_gaussian (image  *g) 

{ 

Pyramid  *p; 

p  =  (Pyramid  * )malloc (sizeof (Pyramid) ) ; 

p->10  =  resize (g,  5) ; 
p->Il  =  reduce (p->10) ? 
p->12  =  reduce (p->Il) ; 
p->13  =  reduce (p->12) ; 
p->l4  =  reduce (p->13 ) ; 
p->15  =  reduce (p->l4) ? 

return(p); 

} 


#endif 
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/**********************************  leplace_n . cpp  ******************************/ 

/*  This  file  contains  the  code  for  generating  a  Leplacian  pyramid  for  a  */ 

/*  given  image.  The  pyramid  depth  is  six  levels.  */ 

/★★♦★★★★A***********************************************************************/ 

#ifndef  LEPLACE 
#define  LEPLACE 

#include  " thesis. h" 

# include  "gaussian__n.  cpp" 

/*  Type  Definitions  */ 

/*  Image  Processing  */ 

Pyramid  *  gen_py r ami d_lep lace  (image  *g) ; 
image  *  cons true t_leplace (Pyramid  *p) ? 

/**********************************************************************************/ 
/*  Leplace  */ 

/★a********************************************************************************/ 
Pyramid  *  gen__pyramid_leplace  ( image  *i) 

{ 

Pyramid  *g,  *p; 
image  *  temp; 

p  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 
g  =  gen  pyramid  gaussian(i) ; 

p->15  =  copy_image (g->15) ; 

temp  =  expand (g->15 ) ; 

p->l4  =  subtract_sign(g->14,  temp) ; 

f ree_image ( temp ) ; 

temp  =  expand (g->l4 ) ; 

p->13  =  subtract_sign (g->l3 ,  temp); 

free_image (temp) ; 

temp  =  expand (g-> 13 ) ; 

p->12  =  subtract_sign (g->12,  temp) ; 

f ree_image ( temp ) ; 

temp  =  expand (g->l2 ) ; 

p->Il  =  subtract_sign (g->Il,  temp) ; 

f ree_image ( temp ) ; 

temp  =  expand (g- >11)  ; 

p->10  =  subtract_sign(g->10/  temp); 

free_image (temp) ; 

return (p) ; 

} 


image  *  c ons t rue t_lep lace (Pyramid  *L) 

{ 

image  *temp,  *i; 

Pyramid  *p; 

p  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

temp  =  expand (L->I5) ;  p->14  =  add_sign (temp,  L->I4);  free (temp); 

temp  =  expand (p->14 ) ;  p->l3  =  add_sign (temp,  L->I3);  free(temp); 

temp  =  expand (p->13 ) ;  p->12  =  add_sign (temp,  L->I2);  free(temp); 

temp  =  expand (p->12 ) ;  p->Il  =  add_sign ( temp,  L->I1) ;  free(temp); 

temp  =  expand (p->Il ) ;  p->10  =  add_sign (temp,  L->I0);  free(temp); 

i  =  copy_image(p->10) ; 

return ( i ) ; 

} 


#endif 
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/***********************************  RoLP„n.cpp  *********************************/ 

/*  This  file  contains  the  code  for  generating  a  RoLP  pyramid  for  a  given  */ 

/*  image.  The  pyramid  depth  is  six  levels.  */ 

^*******************************************************************************/ 
#ifndef  ROLP 
#def ine  ROLP 

#include  "gaussian_n. cpp' 

Pyramid  *  gen_pyramid_RoLP ( image  *g)  ; 
image  *  construct_RoLP (Pyramid  *p) ; 
image  *  max_RoLP { image  *il,  image  *i2); 

/*****************  a****************************************************************/ 

/*  RoLP  */ 

/★★a*******************************************************************************/ 

/*  This  functin  generates  a  six  level  RoLP  pyramid  for  the  given  image  */ 

Pyramid  *  gen_pyramid„RoLP { image  *i) 

{ 

Pyramid  *p,  *g; 
image  *  temp; 

p  =  (Pyramid  *)malloc (sizeof (Pyramid) ) ; 
g  =  gen  _pyramid_gaussian(i) ; 

p->15  =  copy_image (g->15)  ; 

temp  =  expand (g->15 ) ;  p->14  =  divide (g->l4 ,  temp);  f ree_image ( temp ) ; 

temp  =  expand (g->l4 ) ;  p->13  =  divide (g->l3 ,  temp);  f ree_image ( temp ) ; 

temp  =  expand (g->l3 ) ;  p->l2  =  divide (g->l2 ,  temp);  f re e_ image (temp) ; 

temp  =  expand (g->l2 ) ;  p->H  =  divide (g->Il ,  temp);  f ree_ image (temp) ; 

temp  =  expand (g- >11 ) ;  p->lO  =  divide (g->10 ,  temp);  free_image (temp) ; 

return (p) ; 

> 

/*  This  function  returns  the  source  image  from  a  RoLP  pyramid.  */ 
image  *  cons true t_RoLP (Pyramid  *r) 

{ 

image  *  temp ,  * i ; 

Pyramid  *p; 

p  =  (Pyramid  *)malloc (sizeof (Pyramid) ) ; 
i  =  (image  *)malloc (sizeof (image) ) ; 

if (r->15  ! =  NULL)  { 

p->15  =  copy_image (r->15)  ; 
temp  =  expand (p-> 15) ; 
p->14  =  multiply (r->14,  temp); 
free (temp) ; 

} 

else  { 

p->14  =  copy_image(r->I4) ; 

> 

temp  =  expand (p->14) ; 

p->13  =  multiply (r->13 ,  temp); 

free (temp) ; 

temp  =  expand (p->l3 ) ; 

p->12  =  multiply (r->l2 ,  temp); 

free(temp); 

temp  =  expand (p->12 ) ; 

p->Il  =  multiply (r->Il#  temp); 

free (temp) ; 

temp  =  expand (p->Il) ; 

p->lO  =  multiply (r->10 ,  temp) ; 

i  =  multiply (r->10 ,  temp); 

free (temp) ; 

return ( i ) ; 


#endif 
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/it**********************************  gradient_n . cpp  ****************************/ 

/*  This  file  contains  the  code  for  generating  a  gradient  pyramid  for  a  given  */ 

/*  image.  The  pyramid  depth  is  six  levels.  The  gradient  pyramid  is  */ 

/*  actually  four  pyramids  so  a  new  structure  is  also  defined  and  used.  */ 

/***★***************************************★***********************************/ 

#ifndef  GRADIENT 
#def ine  GRADIENT 

#include  " thesis. h" 

# include  ■ convolve_n . cpp " 

# inc lude  " gaus s i an_n . cpp " 

/*  Type  Definitions  */ 

typedef  struct  Pyramid_Gradient  { 

Pyramid  *g,  *dl,  *d2,  *d3,  *d4; 

}  Pyramid_Gradient ; 

/*  Image  Processing  */ 

void  free_Pyramid_Gradient (Pyramid_Gradient  *p)  ; 

Pyramid_Gradient  *  gen_pyramid_gradient (image  *g) ; 
image  *  construct_gradient  (Pyramid__Gradient  *p)  ; 

/*******★★***★**********************★**********************************************/ 

/*  Gradient  */ 

/★it********************************************************************************/ 

/*  This  function  frees  the  memory  used  for  a  gradient  pyramid.  */ 
void  free_Pyramid_Gradient ( Pyramid_Gradient  *p) 

{ 

if (p  ! =  NULL)  { 

free_pyramid (p->g) ; 
free_pyramid{p->dl) ; 
free_pyramid(p->d2) ; 
f ree_pyramid (p->d3 ) ; 
free_pyramid (p->d4) ; 
free (p) ; 
p  =  NULL; 

} 

return; 

> 

/*  This  function  generates  the  four  seperate  pyramids  which  for  the  gradient  pyramid.  */ 
void  gen_gradient (image  *g,  image  **gl,  image  **g2,  image  **g3,  image  **g4) 

{ 

image  *  temp,  *temp2; 

temp  =  conv_w_g_dot„image (g) ; 
temp2  =  add_sign (temp,  g) ; 
f ree_image ( temp ) ; 

( *gl )  =  conv_dl_image ( temp2 ,  1 ) ; 

( *g2 )  =  conv_d2_image ( temp2 ,  1 ) ; 

{ *g3  )  =  conv__d3_image  ( temp2  ,  1 )  ; 

( *g4 )  =  conv_d4_image ( temp2 ,  1 ) ; 
f ree_image ( temp2 ) ; 

return; 

> 

/*  This  function  generates  the  gradient  pyramid  for  a  given  image.  */ 

Pyramid„Gradient  *  gen_py r am id_gradient (image  *i) 

{ 

Pyramid_Gradient  *p; 

p  =  (Pyramid_Gradient  *)malloc (sizeof (Pyramid_Gradient) ) ; 

p->dl  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

p->d2  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

p->d3  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

p->d4  =  (Pyramid  *) malloc (sizeof (Pyramid) ) ; 

p->g  =  (Pyramid  *)malloc (sizeof (Pyramid) ) ; 
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p->g  =  gen_pyramid_gaussian(i) ; 

gen  gradient (p->q-> 10 ,  &p->dl->10,  &p->d2->10,  &p->d3->l0,  &p->d4->l0) 
gen_gradient (p->g->Il ,  &p->dl->Il,  &p->d2->Il,  &p->d3->Il,  &p->d4->Il) 
gen  gradient (p->g->12 ,  &p->dl->12,  &p->d2->12,  &p->d3->l2,  &p->d4->l2) 
gen  gradient (p->q->!3 ,  &p->dl->l3,  &p->d2->l3,  &p->d3->l3,  &p->d4~>13) 
gen  gradient (p->q->14 ,  &p->dl->14,  &p->d2->14,  &p->d3->14,  &p->d4->14) 
gen_gradient (p->g->15 ,  &p->dl->15,  &p->d2->15,  &p->d3->15,  &p->d4->l5) 
return (p) ; 

> 

/*  This  function  returns  the  source  image  from  a  gradient  pyramid.  */ 
image  *  construct_gradient (Pyramid_Gradient  *p) 

{ 

image  *g3 ,  *g2,  *gl,  *g0,  *g4; 
image  *13,  *12,  *11,  *10,  *14; 
image  *1x1,  *1x2,  *1x3,  *1x4; 
image  *lfl,  *lf2,  *lf3,  *lf0,  *lf4; 

lxl  =  gen_gradient_FSD_dl (p->dl->l4) ; 

1x2  =  gen__gradient_FSD„d2  (p->d2->14)  ; 

1x3  =  gen_gradient_FSD_d3 (p->d3->14) ; 

1x4  =  gen_gradient_FSD_d4 (p->d4->l4) ; 
lf4  =  addlltol4 (lxl,  1x2,  1x3,  1x4); 

f ree_image ( lxl ) ;  f ree_image ( 1x2 ) ;  f ree„image ( 1x3 ) ;  f ree_image ( 1x4 ) ? 

lxl  =  gen_gradient_FSD_dl (p->dl->13 ) ; 

1x2  =  gen_gradient_FSD_d2 (p->d2->13 ) ; 

1x3  =  gen_gradient_FSD_d3 (p->d3->13 ) ; 

1x4  =  gen_.gr adient_FSD_d4  (p->d4->13  )  ; 
lf3  =  addllto!4 (lxl,  1x2,  1x3,  1x4); 

f ree_image ( lxl ) ;  f ree_image ( 1x2 ) ;  f ree_image ( 1x3 ) ;  f ree_image ( 1x4 ) ; 

lxl  =  gen_gradient_FSD_dl (p->dl->l2) ; 

1x2  =  gen_gradient_FSD_d2 (p->d2->12 ) ; 

1x3  =  gen_gradient_FSD_d3 (p->d3->12) ; 

1x4  =  gen_gradient_FSD_d4 (p->d4->l2) ; 
lf2  =  addl 1 tol4 { lxl ,  1x2,  1x3,  1x4); 

f ree_image ( lxl ) ;  f ree_image ( 1x2 ) ;  f ree_image ( 1x3 ) ;  f ree_image ( 1x4 ) ; 

lxl  =  gen  gradient  FSD _dl (p->dl->Il) ; 

1x2  =  gen_gradient_FSD_d2 (p->d2->Il) ; 

1x3  =  gen_gradient_FSD_d3 (p->d3->Il) ; 

1x4  =  gen_gradient__FSD_d4  (p->d4->Il)  ; 
lfl  =  addl ltol4 (lxl,  1x2,  1x3,  1x4); 

f ree_image ( lxl ) ;  f ree_image ( 1x2 ) ;  f ree_image ( 1x3 ) ;  f ree_image ( 1x4 ) ; 

lxl  =  gen_gradient_FSD_dl (p->dl->10) ; 

1x2  =  gen_gradient_FSD_d2 (p->d2->l0) ; 

1x3  =  gen  gradient  FSD_.d3  (p->d3->10)  ; 

1x4  =  gen_gradient_FSD_d4 (p->d4->10 ) ; 
lfO  =  addlltol4 ( lxl ,  1x2,  1x3,  1x4); 

f ree_image ( lxl ) ;  f ree_image ( 1x2 ) ;  f ree_image ( 1x3 ) ;  f ree_image ( 1x4 ) ; 

14  =  FSD_to„RE ( If 4 ) ; 

13  =  FSD_to_RE ( If 3 ) ; 

12  =  FSD_to_RE ( If 2 ) ; 

11  =  FSD_tO_RE (lfl) ; 

10  =  FSD_to_RE (lfO ) ; 

g4  =  add_sign(14,  expand_g (p->g->15) ) ; 
g3  =  add_sign ( 13 ,  expand_g (g4 ) ) ; 
g2  =  add_sign ( 12 ,  expand_g ( g3 ) ) ; 
gl  =  add_s ign (11,  expand_g ( g2 ) ) ; 
gO  =  add_sign ( 10 ,  expand_g (gl ) ) ; 

f ree_image (gl ) ;  f ree_image (g2 ) ;  f ree_image (g3 ) ;  f ree_image (g4 ) ; 
f ree_image ( 13 ) ;  f ree_image ( 12 ) ;  f ree_image ( 11 ) ;  f ree_image ( 10 ) ; 
return ( gO ) ; 

} 

#endif 
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/**********************************  m0rph_n.cpp  ************************★★******/ 
/*  This  file  contains  the  code  for  generating  a  morphological  pyramid  for  a  */ 

/*  given  image.  The  pyramid  depth  is  six  levels.  */ 

/***********************★****************************************★**************/ 
#ifndef  MORPH 
#def ine  MORPH 

#include  * thesis. h" 

# include  "files.cpp" 

# include  ■  image_jp .  cpp " 

/*  Type  Definitions  */ 
typedef  struct  Struct_Elem  { 
int  xsize,  ysize,  xcen,  ycen; 

Pixel  *pixel; 

}  Struct_Elem; 

/*  Image  Processing  */ 

Pyramid  *  gen_pyramid_morph ( image  *g)  ; 
image  *  cons t rue t_morph (Pyramid  *p) ; 
image  *  exp and_m ( image  *g) ; 
image  *  reduce_m( image  *g) ; 

/*  Global  Variables  */ 

Struct_Elem  K; 

/*******************************************************************************/ 
/*  STRUCTURING  ELEMENT  STUFF  */ 

/★a*****************************************************************************/ 

/*  Defines  the  structuring  element  as  a  3-by-3  brick.  */ 
void  init_K() 

{ 

K.xsize  as  3;  K.ysize  =  3; 

K.ycen  =  1;  K.xcen  =  1; 

K. pixel  =  (Pixel  *)malloc (K.xsize*K.ysize*sizeof (Pixel) ) ; 


K. pixel [0] [0]  =  1;  K. pixel [1] [0]  =  1;  K. pixel [2] [0]  =  1 

K. pixel [3] [0]  =  1;  K. pixel [4] [0]  =  1;  K. pixel [5] [0 ]  =  1 

K. pixel [6] [0]  =  1;  K. pixel [7] [0]  =  1;  K. pixel [8] [0]  =  1 

K. pixel [0] [1]  =  1;  K. pixel [1] [1]  =  1;  K. pixel [2] [1]  =  1 

K. pixel [3 ] [1]  =  1;  K. pixel [4] [1]  =  1;  K. pixel [5] [1]  =  1 

K. pixel [6] [1]  =  1;  K. pixel [7 ] [1]  =  1;  K. pixel [8] [1]  =  1 

K. pixel [0] [2]  =  1;  K. pixel [1] [2 ]  =  1;  K. pixel [2] [2]  =  1 

K. pixel [3 ] [2]  =  1;  K. pixel [4] [2]  =  1;  K. pixel [5] [2 ]  =  1 

K. pixel [6]  [2  3  =  1;  K. pixel [7]  [2]  =  1;  K. pixel [8]  [2]  =  1 


} 

/*  Defines  the  structuring  element  as  a  2-by-2  brick.  */ 
void  init_ sK() 

{ 

K.xsize  =  2;  K.ysize  =  2; 

K.ycen  =  1;  K.xcen  =  1; 

K. pixel  =  (Pixel  *)malloc(K.xsize*K.ysize*sizeof(Pixel)); 
K. pixel [0] [0 ]  =1;  K. pixel [1) [0]  =  1; 

K. pixel [2] [0]  =  1;  K. pixel [3 ] [0]  =  1; 

K. pixel [0] [1]  =  1;  K. pixel [1] [1]  =  1; 

K. pixel [2] [1]  =  1;  K. pixel [3 ] [1]  =  1; 

K. pixel [0] [2]  =  1;  K. pixel [ 1] [2]  =  1; 

K. pixel [2] [2]  =  1?  K. pixel [3] [2]  =  1; 

} 

/*  Defines  the  structuring  element  as  a  10-by-10  brick.  */ 
void  init_bigK ( ) 

{ 

int  index; 

K.xsize  =  10;  K.ysize  =  10; 


97 


K.ycen  =  5;  K.xcen  =  5; 

K. pixel  =  (Pixel  *)malloc (K.xsize*K.ysize*sizeof (Pixel) ) ; 
for(int  i  =  0;  i  <  K.xsize;  i++)  { 

for(int  j  =  0;  j  <  K.ysize;  j++)  { 

index  =  j*K.xsize  +  i; 

K. pixel [index] [0]  =  1;  K. pixel [index] [ 1]  =  1;  K. pixel [index] [2]  =  1; 

} 

} 

} 

/******************************★****★****★**************************************/ 
/*  IMAGE  PROCESSING  STUFF  */ 

f *******************************************************************************/ 

/*  Returns  a  single  pixel  value  at  the  (i,  j)  location  for  an  image  being 
dilated.  */ 

double  s amp le_dil at e ( image  *  f,  Struct_Elem  *k,  int  i,  int  j,  int  c) 

{ 

int  x,  y,  kx,  ky,  m,  n; 
double  sum,  max  =  -1.0; 

for(m  =  0;  m  <  k->xsize;  m++)  { 

for (n  =0;  n  <  k->ysize;  n++)  { 

kx  =  m  -  k->xcen; 
ky  =  n  -  k->ycen; 
x  =  i  -  kx; 
y  =  j  -  ky; 

if (x  <  0  ||  y  <  0  ||  x  >=  f->xsize  | |  y  >=  f->ysize)  continue; 
sum  =  f->pixel [y*f->xsize  +  x] [c]  +  k->pixel [n*k->xsize  +  m] [c] ; 
if  (sum  >  max)  max  =  sum; 

} 

} 

return  (max)  ; 

} 

/*  Dilates  the  image  f  by  the  structuring  element  K.  */ 
image  *  dilate (image  *f,  Struct_Elem  *k) 

{ 

image  *newpic; 
int  x,  y,  index; 

newpic  =  (image  *)malloc (sizeof (image) ) ; 
newpic->xsize  =  f->xsize; 
newpic->ysize  =  f->ysize; 

newpic->pixel  =  (Pixel  *)malloc (newpic->xsize*newpic->ysize*sizeof (Pixel) ) ; 
for(y  =  0;  y  <  newpic->ysize;  y++)  { 

for(x  =  0;  x  <  newpic->xsize;  x++)  { 

index  =  y*newpic->xsize  +  x; 

newpic->pixel [index] [0]  =  sample_dilate (f ,  k,  x,  y,  0); 

newpic->pixel [index] [1]  =  sample_dilate (f ,  k,  x,  y,  1); 

newpic ->pixel [index] [2]  =  sample_dilate (f ,  k,  x,  y,  2); 

} 

} 

return (newpic) ; 


/*  Returns  a  single  pixel  value  at  the  (i,  j)  location  for  an  image  being 
eroded.  */ 

double  s amp 1 e_erode ( image  *  f,  Struct_J3lem  *k,  int  i,  int  j,  int  c) 

{ 

int  x,  y,  m,  n,  kx,  ky; 
double  sum,  min  =  1000000; 

for(m  =  0;  m  <  k->xsize;  m++)  { 

for (n  =0;  n  <  k->ysize;  n++)  { 

kx  =  m  -  k->xcen ; 
ky  =  n  -  k->ycen; 
x  =  i  +  kx; 
y  =  j  +  ky; 

if (x  <  0  ||  y  <  0  ||  x  >=  f->xsize  | |  y  >=  f->ysize)  continue; 
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sum  =  f->pixel  [y*f->xsize  +  x] [c]  -  k->pixel [n*k~>xsize  +  m]  [c]  ; 
if (sum  <  min)  min  =  sum; 

} 

> 

return (min) ; 

} 

/*  Erodes  the  image  f  by  the  structuring  element  K.  */ 
image  *  erode (image  *f ,  Struct_Elem  *k) 

{ 

image  * newpic; 
int  x,  y,  index; 

newpic  =  (image  *)malloc (sizeof (image) ) ; 
newpic->xsize  =  f->xsize; 
newpic->ysize  =  f->ysize; 

newpic->pixel  =  (Pixel  * )malloc (newpic->xsize*newpic->ysize*sizeof (Pixel) ) ; 
for(y  =  0;  y  <  newpic->ysize;  y++)  { 

for(x  =  0;  x  <  newpic- >xsize;  x++)  { 

index  =  y*newpic->xsize  +  x; 

newpic->pixel [index] [0]  =  sample_erode (f ,  k#  x,  y#  0)  ; 

newpic->pixel [index] [1 3  =  samp le„er ode (f,  k,  x,  y,  1); 

newpic->pixel [index] [2]  =  sample_erode (f ,  k,  x#  y,  2); 

} 

} 

return (newpic) ; 

} 

/*******************************************************************************/ 
/*  FILTERS  */ 

/★★★★★★A************************************************************************/ 

/*  Returns  the  image  f  after  the  application  of  the  opening  morphological 
filter  that  uses  the  structuring  element  K.  */ 
image  *  open (image  *f) 

{ 

image  *temp,  *temp2; 
temp  =  erode  (f,  &K)  ; 
temp2  =  dilate (temp,  &K) ; 
free_image (temp) ; 
return (temp2 ) ? 

} 

/*  Returns  the  image  f  after  the  application  of  the  closing  morphological 
filter  that  uses  the  structuring  element  K.  */ 
image  *  close (image  *f) 

{ 

image  ♦temp,  *temp2; 

temp  =  dilate (f,  &K) ;  //  savef ile ( temp,  *g3rod.ppm’ ) ; 

temp2  =  erode(temp,  &K) ;  //  savef ile (temp2 ,  "g3rode .ppm" ) ; 
f ree_image ( temp) ; 
return ( temp2 ) ; 

> 

/*  Returns  the  image  f  after  the  application  of  both  the  closing  and  opening 
morphological  filters.  */ 
image  *  FCO ( image  *f) 

{ 

image  *temp,  *temp2; 

temp  =  close ( f ) ;  //  savef ile (temp,  wg3ec_b.ppm* ) ; 

temp2  =  open(temp);  //  savef ile (temp2 ,  "g3eo_b.ppm" ) ? 
free_image (temp) ; 
return ( temp2 ) ? 

> 

/*  Returns  the  image  f  after  the  application  of  both  the  opening  and  closing 
morphological  filters.  */ 
image  *  FOC ( image  *f) 

{ 

image  *temp,  *temp2; 

temp  =  open(f);  //  savef ile (temp,  "g3o_b.ppm" ) ; 
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temp2  =  close{temp);  //  savefile (temp2 ,  ’ g3 oc_b . ppm " ) ; 
free_image (temp) ; 
return (temp2) ; 

} 

/*  Returns  the  image  f  after  an  expand  operation  and  closing  filter.  */ 
image  *  C_exp and ( image  *f) 

{ 

image  *temp,  *temp2; 
temp  -  expand_m(f) ; 
temp2  =  close (temp); 
f ree_image ( temp ) ; 
return { temp2 ) ? 

} 

/*  Returns  the  image  f  after  an  expand  operation  and  dilation.  */ 
image  *  D_expand ( image  *f) 

{ 

image  *  temp ,  *  temp2 ; 
temp  =  expand_m ( f ) ; 
temp2  =  dilate (temp,  &K) ; 
f ree_image ( temp ) ; 
return ( temp2 ) ; 

> 

/*  Returns  the  image  f  after  an  expand  operation  and  a  closeing 
then  opening  filter.  */ 
image  *  FCO_expand ( image  *f) 

{ 

image  *temp,  *temp2; 
temp  =  expand_m ( f ) ? 
temp2  =  FCO(temp); 
f ree_image ( temp ) ; 
return (temp2 ) ; 


/*  Returns  the  image  f  after  an  expand  operation  and  an  opening 
then  closeing  filter.  */ 
image  *  FOC_exp and ( image  *f) 

{ 

image  *temp,  *temp2; 

temp  =  FOC { f ) ;  //  savefile (temp,  wg3foc .ppm" ) ; 

temp2  =  exp and_m (temp) ;  //  savefile (temp2 ,  "g3foce.ppm" ) ; 
f ree_image ( temp ) ; 
return ( temp2 ) ; 

} 

/*  Generates  the  next  higher  level  of  a  morphological  pyramid.  */ 
image  *  PC ( image  *  f ) 

{ 

image  *temp,  *temp2; 
temp  =  reduce_m { f ) ? 
temp2  =  FCO(temp); 
f ree_image ( temp ) ; 
return { temp2 ) ; 

} 

/*  Used  for  pyramid  reconstruction  to  generate  the  level  below  the 
current  level.  */ 
image  * PRC (image  *f,  image  *g) 

{ 

image  *temp,  *temp2; 

temp  =  C_expand ( f ) ; 

temp 2  =  subtract_sign (g,  temp); 

f ree_image ( temp ) ; 

return ( temp2 ) ; 
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/*  RESIZING  ROUTINES  */ 

/***************★***************************************************************/ 
/*  Subsamples  the  image.  */ 
image  *  reduce„m ( image  *g) 

{ 

image  * newpic; 

int  xnew,  ynew,  x,  y,  index,  i2; 
newpic  =  (image  *)malloc (sizeof (image) ) ; 
xnew  =  (g->xsize  -  l)/2  +  1; 
ynew  =  (g->ysize  -  l)/2  +  1; 

newpic- >pixel  =  (Pixel  * )malloc (xnew*ynew*sizeof (Pixel) ) ; 
for(x  =  0;  x  <  xnew;  x++)  { 

for (y  =  0;  y  <  ynew;  y++)  { 

index  =  y*xnew  +  x; 
i2  =  y*2*g->xsize  +  x*2; 

newpic->pixel [index] [0]  =  g->pixel[i2] [0]; 
newpic->pixel [index] [1]  =  g->pixel [i2 ] [1] ; 
newpic ->pixel [index] [2]  =  g->pixel [i2] [2] ; 

} 

} 

newpic->xsize  =  xnew; 
newpic ->ysize  =  ynew; 
return (newpic) ; 


/*  Upsamples  the  image.  */ 
image  *  exp and_m ( image  *g) 

{ 

image  *  ne wp i c ; 

int  xnew,  ynew,  x,  y,  index,  i2; 
newpic  =  (image  *)malloc (sizeof (image) ) ; 
xnew  =  (g->xsize  -  1)*2  +  1; 
ynew  =  (g->ysize  -  1)*2  +  1; 

newpic->pixel  =  (Pixel  *)malloc (xnew*ynew* sizeof (Pixel) ) ? 
for(x  =  0;  x  <  xnew;  x++)  { 

for (y  =  0;  y  <  ynew;  y++)  { 

index  =  y*xnew  +  x; 
if (x%2  ==  0  &&  y%2  ==  0)  { 

i2  =  y*g->xsize/2  +  x/2; 

newpic->pixel [index] [0]  =  g->pixel [i2] [0] ; 
newpic->pixel [index] [1]  =  g->pixel[i2][l]; 
newpic->pixel [index] [2]  =  g->pixel [ i2] [2] ; 


} 

else  { 

newpic->pixel [index] [0]  =  0; 
newpic->pixel [index] [1]  =  0; 
newpic->pixel [index] [2]  =  0; 

} 


} 


} 

newpic->xsize  =  xnew; 
newpic ->ysize  =  ynew; 
return (newpic) ? 


/★*********************************************★*******★*************★***★******★**/ 
/*  Morph  */ 

/★★a*********#*********************************************************************/ 

/*  This  functin  generates  a  six  level  morphological  pyramid  for  the  given  image  */ 
Pyramid  *  gen_pyramid_morph ( image  *g) 

{ 

Pyramid  *p; 

image  *  temp,  *i0,  *il,  *i2,  *i3,  *i4,  *i5; 


p  =  (Pyramid  *) malloc (sizeof (Pyramid)); 


//A  regular  resize  instead  of  resize_nn  may  be  appropriate. 

//  resize_nn  keeps  the  edges  of  the  box  from  bluring. 

temp  =  resize_nn (g,  5);  iO  =  close (temp);  free_image (temp) ; 
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il  =  PC(iO) ; 

// 

savef ile (il. 

"gl.ppm") ; 

i2  =  PC  ( i  1 )  7 

// 

savef ile (i2. 

"g2.ppm") ; 

i3  =  PC (i2 ) 7 

// 

savef ile (i3 , 

"g3 .ppm" ) ? 

i4  =  PC (i3 )  ; 

// 

savefile (i4 , 

"g4 .ppm" ) ; 

i5  =  PC (14) ; 

// 

savefile (15, 

"g5.ppm") ; 

p->l5  =  i5 ; 

p->l4  =  PRC (i5 ,  i4) ;  //  savef ile (p->l4/  "14.ppm"); 

p->l3  =  PRC { i4 ,  13 ) ;  //  savef ile (p->13 ,  "13. ppm"); 

p->12  =  PRC ( 13 ,  i2 ) ;  //  savef ile <p->12 ,  "pi2.ppm"); 

p->Il  =  PRC (12 ,  il);  //  savefile (p->Il,  "11. ppm"); 

p->10  =  PRC (11,  iO);  //  savefile (p->10 ,  "lO.ppm"); 

f ree„image ( 14 ) ;  f ree_image (13) ;  f ree_image ( 12 ) ; 
f ree_image ( i 1 ) ;  f ree_image ( i 0 ) ; 

return (p) ; 

> 

/*  This  function  returns  the  source  image  from  a  morphological  pyramid.  */ 
image  *  c onst rue t_morph (Pyramid  *p) 

{ 

image  *g4e,  *g3e,  *g2e,  *gle,  *g5e; 
image  *i0,  *il,  *i2,  *i3,  *i4; 

g5e  =  C_expand (p->15) ;  //  savefile (g4e,  "g4ec.ppm"); 

i4  =  add_sign (g5e,  p->l4) ;  //  savefile (i3,  "i3c.ppm"); 

g4e  =  C_expand(i4) ;  //  savef ile (g4e,  "g4ec.ppm"); 

i3  =  add_sign (g4e,  p->l3);  //  savefile (i3,  "i3c.ppm"); 

g3e  =  C_expand(i3 ) ;  //  savef ile (g3e,  "g3ec.ppm"); 

i2  =  add_s ign ( g3 e ,  p->l2);  //  savefile(i2,  "i2n.ppm"); 

g2e  =  C_expand ( i2 ) ;  //  savef ile (g2e,  "g2en.ppm"); 

il  =  add_sign (g2e#  p->Il) ;  //  savefile(il,  "iln.ppm"); 

gle  =  C_expand(il) ?  //  savef ile (gle,  "glen. ppm"); 

iO  =  add_sign(gle/  p->10);  //  savefile (iO,  "iOn.ppm"); 

free_image (gle) ;  free_image (g2e) ;  free_image (g3e) ;  free_image (g4e) ; 

f ree_image ( 13 ) ;  f ree_image ( i2 ) ;  f ree_image ( il ) ; 

return (iO) ; 

} 


#endif 
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/*******★*************************  convolve. cpp  ********************************/ 
/*  This  file  contains  the  code  for  implementing  convolution  for  pyramid  */ 

/*  filtering.  It  also  contains  code  necessary  to  support  the  gradient  */ 

/*  pyramid  */ 

/♦★★a***************************************************************************/ 

#ifndef  CONVOLVE 
#def ine  CONVOLVE 

#include  ’ thesis. h" 

#include  'files. cpp* 

# inc lude  " image_p . cpp  * 

#def ine  sqrt2  1.41421356237 
#def ine  sqrt2inv  .707106781188 
#define  eighth  .125 

/*  Image  Processing  */ 
image  *  reduce_g (image  *g)  ; 
image  *  expand  cr( image  *g); 

/★a*****************************************************************************/ 

/*  GAUSSIAN  KERNEL  FUNCTIONS  */ 

/****★**************************************************************************/ 
/*  Gaussian  kernal  */ 
double  w_g(int  m,  int  n) 

{ 

double  wm,  wn,  answer; 
if (m  ==  0)  wm  =  6; 

else  if  (abs(m)  ==  1)  wm  =  4; 

else  wm  =  1; 

if (n  ==  0 )  wn  =  6; 

else  if  (abs(n)  ==  1)  wn  =  4; 

else  wn  =  1; 

answer  =  wm*wn/256; 
return (answer) ; 

} 

/*  Seperated  function  of  Gaussian  kernal  */ 
double  w_g_dot(int  m#  int  n) 

{ 

double  wm,  wn,  answer; 


if (m  == 

0) 

wm  =  2 

else 

wm  =  1 

if (n  == 

0) 

wn  =  2 

else 

wn  =  1 

answer  =  wm*wn/16; 
return (answer) ; 

} 

/*********★****************************************************★******★*********/ 
/*  CONVOLUTION  ROUTINES  */ 

/★★★★★★★A***********************************************************************/ 

/*  Returns  a  single  sample  located  at  (i,  j)  for  a  given  color  (red,  green,  or 
blue  pixels)  for  the  convolution  of  image  with  w_dot.  */ 
double  sample  w  g  dot ( image  *  pic,  int  i,  int  j,  int  color) 

{ 

int  m,  n,  index,  tx,  ty; 
double  tsum  =  0.0; 
for(m  =  -1;  m  <  2;  m++)  { 

for(n  =  -1;  n  <  2;  n++)  { 

tx  =  i  +  m; 
if(tx  <  0)  tx  =  -tx; 

if(tx  >=  pic->xsize)  tx  =  2*pic->xsize  -  tx  -  1; 
ty  =  j  +  n; 


} 


if (ty  <  0)  ty  =  -ty; 

if(ty  >=  pic->ysize)  ty  =  2*pic->ysize  -  ty  -  1; 
index  =  ty*pic->xsize  +  tx; 

tsum  +=  w_g_dot(m,n)  *  pic->pixel [index] [color] ; 


> 

return  (tsum)  ; 


/*  Returns  a  single  sample  located  at  (i,  j )  for  a  given  color  (red,  green,  or 
blue  pixels)  for  the  convolution  of  image  with  w  the  gaussian  kernel.  */ 
double  s amp le_w_g (image  *  pic,  int  i,  int  j ,  int  color) 


int  m,  n,  index,  tx,  ty; 
double  tsum  =  0.0; 
for(m  =  -2;  m  <  3;  m++)  { 

for(n  =  -2;  n  <  3;  n++)  { 

tx  =  i  +  m; 
if(tx  <  0)  tx  =  -tx; 

if(tx  >=  pic->xsize)  tx  =  2*pic->xsize  -  tx  -  1; 


ty  =  j  +  n; 

if ( ty  <  0 )  ty  =  -ty; 

if (ty  >=  pic->ysize)  ty  =  2*pic->ysize  -  ty  -  1; 


index  =  ty*pic->xsize  +  tx; 

tsum  +=  w_g(m,n)  *  pic->pixel [index] [color] ; 

} 

} 

return (tsum)  ; 

} 

/*  Returns  a  single  sample  located  at  (i,  j)  for  a  given  color  (red, 
green,  or  blue  pixels)  for  the  convolution  of  image  with  oriented 
derivative  filter  dl .  */ 

double  s amp le_dl (image  *  pic,  int  i,  int  j,  int  color,  double  mult) 

{ 

int  index; 


if(i  >=  0  &&  i  <=  pic->xsize)  { 
index  =  j *pic->xsize  +  i; 

if(i  ==  0)  return (pic->pixel [index] [color] ) *mult/3 . 5; 

else  if(i  ==  pic->xsize)  return (- (pic->pixel [index- 1] [color] ) *mult/3 . 5) ; 
else  return ( (pic->pixel [index] [color]  -  pic->pixel [index- 1] [color] ) *mult) ; 

} 

else 

return (-111000) ; 


/*  Returns  a  single  sample  located  at  (i,  j)  for  a  given  color  (red, 
green,  or  blue  pixels)  for  the  convolution  of  image  with  oriented 
derivative  filter  d2 .  */ 

double  s amp le_d2 (image  *  pic,  int  i,  int  j,  int  color,  double  mult) 

{ 

double  a,  b; 

if (i  >=  0  &&  i  <=  pic->xsize)  { 

if(i  >  0  &&  i  <=  pic->xsize  &&  j  <  pic->ysize)  { 
a  =  pic->pixel [ j *pic->xsize  +  i  -  1] [color] ; 

> 

else  a  =  0; 

if ( j  >  0  &&  j  <=  pic->ysize  &&  i  <  pic->xsize)  { 
b  =  pic->pixel [ ( j-1) *pic->xsize  +  i] [color] ; 

> 

else  b  =  0; 

if (a  ==  0  | |  b  ==  0) 

return((b  -  a) *sqrt2inv*mult/3 . 5) ; 
else 

return ( (b  -  a) *sqrt2inv*mult) ; 
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} 

else 

return (-111000) ; 

} 

/*  Returns  a  single  sample  located  at  (i,  j)  for  a  given  color  (red, 
green,  or  blue  pixels)  for  the  convolution  of  image  with  oriented 
derivative  filter  d3 .  */ 

double  s amp le_d3 (image  *  pic,  int  i,  int  j,  int  color,  double  mult) 

{ 

double  a,  b; 

if  ( j  >=  0  ScSc  j  <=  pic->ysize)  { 
if ( j  <  pic->ysize) 

a  =  pic->pixel [ j *pic->xsize  +  i] [color] ; 
else  a  =  0; 

if ( j  >  0) 

b  =  pic->pixel [ ( j-1) *pic->xsize  +  i] [color] ; 
else  b  =  0; 
if (a  ==  0  | |  b  ==  0) 

return ( (b  -  a) *mult/3 . 5) ; 
else 

return ( (b  -  a) *mult) ; 

> 

else 

return (-111000 ) ; 


/*  Returns  a  single  sample  located  at  (i,  j)  for  a  given  color  (red, 
green,  or  blue  pixels)  for  the  convolution  of  image  with  oriented 
derivative  filter  d4.  */ 

double  s amp le_d4 (image  *  pic,  int  i,  int  j,  int  color,  double  mult) 

{ 

double  a,  b; 

if(j  >=  0  &&  j  <=  pic->xsize)  { 

if(j  <  pic->ysize  &&  i  <  pic->xsize)  { 
a  =  pic->pixel [ j *pic->xsize  +  i] [color] ; 

> 

else  a  =  0; 

if ( j  >  0  &&  i  >  0)  { 

b  =  pic->pixel [ ( j-1) *pic->xsize  +  (i  -  1) ] [color] ; 

} 

else  b  =  0; 

if (a  ==  0  | |  b  ==  0) 

return ( (b  -  a) *sqrt2inv*mult/3 . 5) ; 
else 

return ( (b  -  a) *sqrt2inv*mult) ; 

} 

else 

return ( -111000 ) ; 

} 

/*  Convolves  the  image  g  with  the  Gaussian  kernel  w  */ 
image  *  conv_w_g_image ( image  *g ) 

{ 

image  *  newp i c ; 

int  xnew,  ynew,  x,  y,  index; 

newpic  =  (image  * )malloc (sizeof ( image) ) ; 

xnew  =  g->xsize; 

ynew  =  g->ysize; 

newpic->pixel  =  (Pixel  *) malloc (xnew*ynew* sizeof (Pixel) ) ; 
for(x  =  0;  x  <  xnew;  x++)  { 

for (y  =  0;  y  <  ynew;  y++)  { 

index  =  y*xnew  +  x; 

newpic ->pixel [index] [0]  =  sample_w_g (g,  x,  y,  0); 
newpic->pixel [index] [1]  =  sample_w_g (g,  x,  y,  1); 
newpic->pixel [index] [2]  =  sample_w„g (g,  x,  y,  2); 

} 


} 

newpic->xsize  =  xnew; 
newpic->ysize  =  ynew; 
return (newpic) ; 

} 

/*  Convolves  the  image  g  with  the  seperated  Gaussian  kernel  w_dot  *7 
image  *  conv_w_g_dot_image ( image  *g) 

{ 

image  *  newp i c ; 

int  xnew,  ynew,  x,  y,  index; 

newpic  =  (image  *) malloc (sizeof (image) ) ; 

xnew  =  g->xsize; 

ynew  =  g->ysize; 

newpic->pixel  =  (Pixel  *) malloc (xnew*ynew*sizeof (Pixel) ) ; 
for(y  =0;  y  <  ynew;  y++)  { 

for(x  =  0;  x  <  xnew;  x++)  { 

index  =  y*xnew  +  x; 

newpic- >pixel [index] [0]  =  sample„w_g_dot (g,  x,  y,  0); 
newpic->pixel [index] [1]  =  sample_w_g_dot (g,  x,  y,  1); 
newpic->pixel [index] [2]  =  sample_w_g_dot (g,  x,  y,  2); 

} 

} 

newpic->xsize  =  xnew; 
newpic->ysize  =  ynew; 
return (newpic) ; 

} 

/*  Convolves  the  image  g  with  the  first  derivative  filter  dl  timesed 
by  a  multiplier  mult  */ 
image  *  conv_dl_image (image  *g,  double  mult) 

{ 

image  *  newpic ; 
int  x,  y,  index; 

newpic  =  (image  *) malloc (sizeof (image) ) ; 
newpic->xsize  =  g->xsize+l; 
newpic->ysize  =  g->ysize; 

newpic ->pixel  =  (Pixel  * )malloc (newpic->xsize*newpic->ysize*sizeof (Pixel ) ) 
for(y  =  0;  y  <  newpic->ysize;  y++)  { 

for (x  =  0;  x  <  newpic->xsize;  x++)  { 

index  =  y*newpic->xsize  +  x; 

newpic->pixel [index] [0 ]  =  sample„dl (g,  x,  y,  0,  mult); 

newpic->pixel [index] [1]  =  sample_dl(g,  x,  y,  1,  mult); 

newpic->pixel [index] [2]  =  sample_dl(g,  x,  y,  2,  mult); 

} 

} 

return (newpic) ; 

} 

/*  Convolves  the  image  g  with  the  first  derivative  filter  d2  timesed 
by  a  multiplier  mult  */ 
image  *  conv_d2_image ( image  *g,  double  mult) 

{ 

image  *  newp i c ; 
int  x,  y,  index; 

newpic  =  (image  *) malloc (sizeof (image) ) ; 
newpic->xsize  =  g->xsize+l; 
newpic~>ysize  =  g->ysize+l; 

newpic->pixel  =  (Pixel  * )malloc (newpic->xsize*newpic->ysize*sizeof (Pixel ) ) 
for(y  =  0;  y  <  newpic->ysize;  y++)  { 

for(x  =  0;  x  <  newpic->xsize;  x++)  { 

index  =  y*newpic->xsize  +  x; 

newpic- >pixel [index] [0 ]  =  sample_d2 (g,  x,  y,  0,  mult); 

newpic ->pixel [index] [1]  =  sample_d2(g,  x,  y,  1,  mult); 

newpic->pixel [index] [2]  =  sample_d2(g,  x,  y,  2,  mult); 

} 

} 

return ( newpic ) ; 
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/*  Convolves  the  image  g  with  the  first  derivative  filter  d3  timesed 
by  a  multiplier  mult  */ 
image  *  conv_d3„image ( image  *g,  double  mult) 

{ 

image  *  newp i c ; 
int  x,  y,  index; 

newpic  =  (image  *) malloc (sizeof ( image) )  ; 
newpic->xsize  =  g->xsize; 
newpic->ysize  =  g->ysize+l; 

newpic->pixel  =  (Pixel  *) malloc (newpic->xsize*newpic->ysize*sizeof (Pixel) ) ; 
for(y  =  0;  y  <  newpic->ysize;  y++)  { 

for(x  =  0;  x  <  newpic->xsize;  x++)  { 

index  =  y*newpic->xsize  +  x; 

newpic- >pixel [index] [0]  =  sample_d3(g,  x,  y,  0,  mult); 

newpic->pixel [index] [1]  =  sample__d3  (g,  x,  y,  1,  mult); 

newpic->pixel [index] [2]  =  sample„d3 (g,  x,  y,  2,  mult); 

} 

} 

return (newpic); 


/*  Convolves  the  image  g  with  the  first  derivative  filter  d4  timesed 
by  a  multiplier  mult  */ 
image  *  conv_d4_image (image  *g,  double  mult) 

{ 

image  * newpic; 
int  x,  y,  index; 

newpic  =  (image  *) malloc (sizeof (image) ) ; 
newpic->xsize  =  g->xsize+l; 
newpic->ysize  =  g->ysize+l; 

newpic->pixel  =  (Pixel  *) malloc (newpic->xsize*newpic->ysize*sizeof (Pixel) ) ; 
for(y  =  0;  y  <  newpic->ysize;  y++)  { 

for(x  =  0;  x  <  newpic->xsize;  x++)  { 

index  =  y*newpic->xsize  +  x; 

newpic->pixel [index] [0]  =  sample_d4(g,  x,  y,  0,  mult); 

newpic->pixel [index] [1]  =  sample_d4(g,  x,  y,  1,  mult); 

newpic- >pixel [index] [2]  =  sample_d4(g,  x,  y,  2,  mult); 

} 

} 

return (newpic) ; 


/********************************************************************** *********/ 
/*  RESIZING  ROUTINES  */ 

/*******************************************************************************/ 
/*  Returns  a  single  sample  located  at  (i,  j)  for  the  REDUCED  image  of  pic 
color  is  the  indicator  for  red,  green,  or  blue  pixels.  */ 
double  reducepixel_g ( image  *  pic,  int  i,  int  j,  int  color) 

{ 

int  m,  n,  index,  tx,  ty; 

double  tsum  =  0.0; 

for(m  =  -2;  m  <  3;  m++)  { 

for(n  =  -2;  n  <  3;  n++)  { 

tx  =  2*i  +  m; 

if(tx  <0  ||  tx  >=  pic->xsize)  continue; 
ty  =  2*j  +  n; 

if(ty  <0  ||  ty  >=  pic->ysize)  continue; 
index  =  ty*pic->xsize  +  tx; 

tsum  +=  w_g(m,n)  *  pic->pixel [index] [color] ; 

> 

} 

return (tsum) ; 


/*  Filters  and  subsamples  the  image  g  to  get  the  next  level  in  the 
gradient  pyramid  */ 


image  *  reduce_g (image  *g) 

{ 

image  *  newp i c ; 

int  xnew,  ynew,  x,  y,  index; 

newpic  =  (image  *) malloc (sizeof (image) ) ; 

xnew  =  (g->xsize  -  l)/2  +  1; 

ynew  =  (g->ysize  -  1 ) / 2  +  1; 

newpic->pixel  =  (Pixel  * )malloc (xnew*ynew*sizeof (Pixel) ) ; 
for(x  =  0;  x  <  xnew;  x++)  { 

for (y  =  0;  y  <  ynew;  y++)  { 

index  =  y*xnew  +  x; 

newpic->pixel [index] [0]  =  (int) reducepixel_g (g,  x,  y,  0); 
newpic->pixel [index] [1]  =  (int) reducepixel_g (g,  x,  y,  1); 
newpic->pixel [index] [2]  =  (int) reducepixel_g (g,  x,  y,  2); 

} 

} 

newpic->xsize  =  xnew; 
newpic->ysize  =  ynew; 
return (newpic) ; 

} 

/*  Returns  a  single  sample  located  at  (i,  j)  for  the  EXPANDED  image  of  pic 
color  is  the  indicator  for  red,  green,  or  blue  pixels.  */ 
double  expandpixel_g (image  *  pic,  int  i,  int  j,  int  color) 

{ 

int  m,  n,  index,  tx,  ty; 
double  tsum  =  0.0; 
for(m  =  -2;  m  <  3;  m++)  { 

for(n  =  -2;  n  <  3;  n++)  { 

if (  ( ( i-m) %2  ==  0)  &&  ( ( j-n) %2  ==  0)  )  { 

tx  =  (i-m)/2; 

if (tx  <0  ||  tx  >=  pic->xsize)  continue; 
ty  =  ( j  -  n)  /2 ; 

if(ty  <0  ||  ty  >=  pic->ysize)  continue; 
index  =  ty*pic->xsize  +  tx; 

tsum  +=  w_g(m,n)  *  pic->pixel [index] [color]; 

} 

} 

} 

return (4 *t sum) ; 


/*  Upsamples  and  interpolates  the  image  g  to  get  a  lower  resolution  image 
of  the  previoius  level  in  the  pyramid*/ 
image  *  expand  q( image  *g) 

{ 

image  *newpic; 

int  xnew,  ynew,  x,  y,  index; 

newpic  =  (image  *) malloc (sizeof (image) ) ; 

xnew  =  (g->xsize  -  1)*2  +  1; 

ynew  =  (g->ysize  -  1)*2  +  1; 

newpic->pixel  =  (Pixel  *) malloc (xnew*ynew*sizeof (Pixel )) ; 
for(x  =  0;  x  <  xnew;  x++)  { 

for(y  =  0;  y  <  ynew;  y++)  { 

index  =  y*xnew  +  x; 

newpic->pixel [index] [0]  =  (int) expandpixel_g (g,  x,  y,  0); 
newpic->pixel [index] [1]  =  (int) expandpixel_g(g,  x,  y,  1); 
newpic->pixel [index] [2]  =  (int) expandpixel_g (g,  x,  y,  2); 

> 

> 

newpic->xsize  =  xnew; 
newpic->ysize  =  ynew; 
return (newpic) ; 

} 

/*  Adds  the  images  of  the  four  different  gradient  pyramids  associated 
with  dl,  d2,  d3,  and  d4.  It  returns  the  FSD  leplacian  pyramid.  */ 
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image  *addlltol4 (image  *11,  image  *12,  image  *13,  image  *14) 

{ 

image  *newpic; 

int  xnew,  ynew,  x,  y,  index; 
int  xtemp,  ytemp; 
xnew  =  13->xsize; 
ynew  =  ll->ysize; 

newpic  =  (image  *) malloc (sizeof (image) ) ; 

newpic->pixel  =  (Pixel  * )malloc (xnew*ynew* sizeof (Pixel) ) ; 
xtemp  =  xnew  +  2; 
for (y  =  0;  y  <  ynew?  y++)  { 

for(x  =  0;  x  <  xnew;  x++)  { 

index  =  y*xnew  +  x; 
ytemp  =  y  +  1? 

newpic- >pixel [index] [0]  =  ll->pixel [y  *xtemp  +  x  +  1]  [0]  + 

12 - >pixel [ytemp *xtemp  +  x  +  1] [0]  + 

13- >pixel [ytemp* (xnew)  +  x  ] [0]  + 

14- >pixel [ytemp*xtemp  +  x  +  1 ] [ 0 ] ; 
newpic- >pixel [index] [1]  =  ll->pixel[y  *xtemp  +  x  +  1] [1]  + 

12- >pixel [ytemp*xtemp  +  x  +  1] [1]  + 

13 - >pixel [ytemp* (xnew)  +  x  ] [1]  + 

14 - >pixel [ytemp *xtemp  +  x  +  1] [1] ; 

newpic ->pixel [index] [2 ]  =  ll->pixel [y  *xtemp  +  x  +  1] [2]  + 

12 - >pixel [ytemp *xtemp  +  x  +  1] [2]  + 

13 - >pixel [ytemp* (xnew)  +  x  ] [2]  + 

14- >pixel [ytemp*xtemp  +  x  +  1 ] [ 2 ] ; 

} 

} 

newpic ->xsize  =  xnew? 
newpic->ysize  =  ynew; 
return (newpic) ? 

> 

/★if**************************************************************************-******/ 

/*  FSD  Conversions  */ 

/***★****************************************************★***************★*********/ 
/*  Converts  the  dl  gradient  image  into  the  dl  FSD  Leplacian  image.  */ 
image  *gen__gradient_FSD_dl  ( image  *d) 

{ 

image  *  temp; 
int  x,  y,  index; 

temp  =  (image  *)malloc (sizeof (image) ) ; 
temp->xsize  =  d->xsize+l; 
temp->ysize  =  d->ysize; 

temp->pixel  =  (Pixel  *) malloc (temp->xsize*temp->ysize*sizeof (Pixel) ) ; 
for(x  =  0;  x  <  temp->xsize;  x++)  { 

for(y  =  0;  y  <  temp->ysize;  y++)  { 

index  =  y*temp->xsize  +  x; 

temp->pixel [index] [0]  =  sample_dl(d,  x,  y,  0,  -eighth); 

temp->pixel [index] [1]  =  sample„dl(d,  x,  y,  1,  -eighth); 

temp->pixel [index] [2]  =  sample_dl (d,  x,  y,  2,  -eighth); 

} 

} 

return ( temp ) ; 


/*  Converts  the  d2  gradient  image  into  the  d2  FSD  Leplacian  image.  */ 
image  *gen  gradient  FSD  d2( image  *g) 

{ 

image  *  temp; 
int  x,  y,  index; 

temp  =  (image  * ) malloc (sizeof ( image) ) ; 
temp->xsize  =  g->xsize  +  1; 
temp->ysize  =  g->ysize  +  1? 

temp->pixel  =  (Pixel  *) malloc (temp->xsize*temp->ysize*sizeof (Pixel )) ; 
for(x  =  0;  x  <  temp->xsize;  x++)  { 

for(y  =  0;  y  <  temp->ysize;  y++)  { 


index  =  y*temp->xsize  +  x; 

temp->pixel [index] [0]  =  sample_d2 (g,  x,  y,  0,  -eighth); 

temp->pixel [index] [1]  =  sample_d2 (g,  x,  y,  1,  -eighth); 

temp->pixel [index] [2 ]  =  sample_d2(g,  x,  y,  2,  -eighth); 

} 

} 

return (temp) ; 

} 

/*  Converts  the  d3  gradient  image  into  the  d3  FSD  Leplacian  image.  */ 
image  *gen_gradient_FSD_d3 { image  *g) 

{ 

image  *  temp; 
int  x,  y,  index; 

temp  =  {image  *)malloc (sizeof (image) ) ; 
temp->xsize  =  g->xsize; 
temp->ysize  =  g->ysize  +  1; 

temp->pixel  =  (Pixel  *) malloc ( temp->xsize*temp->ysize*sizeof (Pixel ) ) 
for (x  =  0;  x  <  temp->xsize;  x++)  { 

for(y  =  0;  y  <  temp->ysize;  y++)  { 

index  =  y*temp->xsize  +  x; 

temp- >pixel [index] [0]  =  sample_d3 (g,  xf  y,  0,  -eighth); 

temp- >pixel [index] [1]  =  sample_d3 (g,  x,  y#  1,  -eighth); 

temp- >pixel [index] [2]  =  sample_d3 (g,  x,  y,  2,  -eighth); 

} 

} 

return (temp) ; 


/*  Converts  the  d4  gradient  image  into  the  d4  FSD  Leplacian  image.  */ 
image  *gen_gradient_FSD_d4 ( image  *g) 

{ 

image  *  temp ; 
int  x,  y,  index; 

temp  =  (image  *) malloc (sizeof (image) ) ; 
temp->xsize  =  g->xsize  +  1; 
temp->ysize  =  g->ysize  +  1; 

temp->pixel  =  (Pixel  *)malloc (temp->xsize*temp->ysize*sizeof (Pixel) ) 
for(x  =  0;  x  <  temp->xsize;  x++)  { 

forty  =  0;  y  <  temp->ysize;  y++)  { 

index  =  y*temp->xsize  +  x; 

temp->pixel [index] [0]  =  sample_d4 (g#  x,  y,  0#  -eighth); 

temp- >pixel [index] [1]  =  sample_d4 (g,  x,  y,  1,  -eighth); 

temp- >pixel [index] [2]  =  sample_d4 (g,  x,  y,  2,  -eighth); 

} 

} 

return (temp) ; 


/*  Converts  the  FSD  Leplacian  image  into  a  Leplacian  image.  */ 
image  *  FSD_to_RE ( image  *g) 

{ 

image  *temp#  *t2; 

t2  =  conv_w„g_image (g) ; 
temp  =  add_sign(g/  t2)  ; 
f ree_image ( t2 )  ; 

return ( temp ) ; 

} 


#endif 


/**********************************  files. cpp  **********************************/ 

/*  This  file  contains  the  code  for  reading  an  image  file,  and  storing  the  */ 

/*  pixel  information  into  an  image  structure.  The  loadfile  routine  is  set  */ 

/*  for  portable  bitmaps,  but  can  be  implemented  for  any  image  storage  type.  */ 

/*  Image  memory  routines  are  also  set  up,  as  well  as  routines  to  resize  the  */ 

/*  image  and  copy  the  image.  (used  for  portable  bitmaps  only)  */ 

/**★****************************************************************************/ 
#ifndef  FILES_CPP 
#def ine  FILES_CPP 

# include  "thesis.h" 

/*  This  is  the  type  of  a  pixel  0  is  red,  1  is  green,  and  index  2  is  blue  */ 
typedef  double  Pixel [3]; 

/*  This  is  the  structure  to  store  an  image  as  a  two  dimentional  array  of  pixels  */ 
typedef  struct  image  { 
int  xsize,  ysize; 

Pixel  *pixel ; 

>  image ; 

/*  Functions  which  are  defined  later  */ 
void  free_image ( image  *  il) ; 

void  initfile (FILE  **infile,  char  *  s,  int  *xs,  int  *ys)  ; 
void  loadfile (image  *filestruct,  char  *  name); 

void  interpolate__pixel_nn (image  *  g,  int  xnew,  int  ynew,  int  x,  int  y,  DPoint  color); 
void  interpolate__pixel_bi  (image  *  g,  int  xnew,  int  ynew,  int  x,  int  y,  DPoint  color); 
image  *  resize (image  *il,  int  n) ; 
image  *  copy_image (image  *i) ; 

/*  Function  definitions  */ 

/*  this  function  copies  the  image  passed  in  to  another  image  */ 

/*  basically  it  allocates  memory  and  copies  the  contents  of  the  image  i  and 
returns  a  pointer  to  the  new  image  */ 
image  *  copy„image ( image  *i) 

{ 

int  index,  x,  y; 
image  *d; 

d  =  (image  *)malloc (sizeof (image) ) ; 

d->pixel  =  (Pixel  * )malloc (i->xsize*i->ysize*sizeof (Pixel )) ; 

d->xsize  =  i->xsize; 
d->ysize  =  i->ysize; 

for(y  =  0;  y  <  i->ysize;  y++)  { 

for(x  =  0;  x  <  i-*>xsize;  x++)  { 

index  =  y*i->xsize  +  x; 

d->pixel [index] [0]  =  (double)i->pixel[index][0]; 
d->pixel [index] [1]  =  (double) i->pixel [index] [1] ; 
d->pixel [index] [2]  =  (double) i->pixel [index] [2] ; 

> 

} 

return (d) ; 

} 

/*  This  frees  the  memory  used  to  store  the  image  */ 
void  free„ image ( image  *  il) 

{ 

if (il  !=  NULL)  { 

if (il->pixel  ! =  NULL)  { 
free (il->pixel) ; 

} 

free (il) ; 

} 

il  =  NULL; 
return; 

> 

/*  reads  the  header  of  a  ppm  file  and  returns  the  x  and  y  dimentions  in  xs  and  ys  */ 
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void  initf ile (FILE  **infile,  char  *  s,  int  *xs,  int  *ys) 

{ 

unsigned  char  dummystring [256] ; 
unsigned  char  final; 
int  xsize,  ysize,  dummy; 

(*infile)  =  fopen(s,  "r")? 
if  ( ! ( *inf ile) ) 

return; 

fscanf ( ( *inf ile) ,  "%s" , dummystring) ; 
fscanf {( *inf ile) ,  "%c" , &f inal ) ; 

while (final  !=  *\n*) 
fscanf ( (*inf ile) ,  *%c* , &final) ; 
fscanf ( (*inf ile) ,  "%c * , &f inal) ; 

while (final  !=  '  \n*) 
fscanf (( *inf ile)  ,  "%c*,&final); 

fscanf ( (*inf ile) ,  "%d%d%d%c" ,  fcxsize,  &ysize,  &dummy,  fcfinal); 
while (final  !=  ’\n') 
fscanf ( (*infile) ,  "%c* , &final) ; 

(*xs)  =  xsize; 

( *ys )  =  ysize; 


/*  this  function  loads  a  file  into  the  image  pointer  filestruct 
name  is  a  string  which  is  the  name  of  the  file  to  be  read  in. 
the  image  must  be  saved  as  a  ppm  to  be  read  using  this  function  */ 
void  loadf ile (image  *filestruct/  char  *  name) 

{ 

int  xsize,  ysize,  x,  y,  index; 

FILE  *  infile; 

int  red,  green,  blue; 

initf ile (&inf ile,  name,  &xsize,  &ysize) ; 

filestruct->pixel  =  (Pixel  *)malloc (xsize*ysize*sizeof (Pixel) ) ; 

for(y  =  0;  y  <  ysize;  y++)  { 

for(x  =  0;  x  <  xsize;  x++)  { 

fscanf (inf ile,  "%d%d%dw,  &red,  &green,  &blue) ; 
index  =  y*xsize  +  x; 
filestruct->pixel [index] [0]  =  red; 
filestruct->pixel [index] [1]  =  green; 
f ilestruct->pixel [index] [2 ]  =  blue; 

} 

} 

f ilestruct->xsize  =  xsize; 

f ilestruct->ysize  =  ysize; 

fclose (infile) ; 


/*  This  function  copies  the  contents  of  the  image  to  a  ppm  file  with  the  name  of 
the  string  passed  in  */ 
void  savef ile (image  *filestruct,  char  *  name) 

{ 

int  numcol  =  255,  x,  y,  xsize,  ysize,  index; 

FILE  *  outfile; 
outfile  =  fopen (name, *ww ) ; 
if  (outfile  ==  NULL) 
return; 

xsize  =  f ilestruct->xsize; 
ysize  =  filestruct->ysize; 

fprintf (outfile,  "P3\n#  Created  by  Ted  Meek\n"); 
fprintf (outfile,  "%d  %d  %d\n",  xsize,  ysize,  numcol); 
for(y  =  0;  y  <  ysize;  y++)  { 

for(x  =  0;  x  <  xsize;  x++)  { 

index  =  y*xsize  +  x; 

fprintf (outfile,  "%d  %d  %d  ",  (int) abs ( f ilestruct~>pixel [index] [0] ) , 

( int) abs ( f ilestruct->pixel [index] [1] ) , 
(int) abs ( f ilestruct->pixel [index]  [2] ) ) ; 


} 


} 

fclose (outfile) ; 


/*  This  function  copies  the  contents  of  the  image  to  a  ppm  file  with  the  name  of 
the  string  passed  in.  It  is  a  grayscale  of  either  the  red,  blue,  or  green  image 
To  save  red,  color  =  0,  green  color  =  1,  and  blue  color  =  2*1 
void  savef ile_onecolor ( image  *filestruct,  char  *  name,  int  color) 

{ 

int  numcol  =  255,  x,  y,  xsize,  ysize,  index; 
int  red,  green,  blue; 

FILE  *  outfile; 
outfile  =  f open (name, "w" ) ; 
if  (outfile  ==  NULL) 
return; 

xsize  =  filestruct->xsize; 
ysize  =  filestruct->ysize; 

f print f (outfile,  "P3\n#  Created  by  Ted  Meek\n»); 
fprintf (outfile,  *%d  %d  %d\n" ,  xsize,  ysize,  numcol); 
for (y  =  0;  y  <  ysize;  y++)  { 

for (x  =  0;  x  <  xsize;  x++)  { 

index  =  y*xsize  +  x; 

red  =  (int ) filestruct->pixel [index] [color] ; 
green  =  (int) f ilestruct->pixel [index] [color] ; 
blue  =  (int) filestruct->pixel [index] [color] ; 
fprintf (outfile,  "%d  %d  %d  ",  red,  green,  blue); 

} 

> 

fclose (outfile) ; 


/*  This  resizes  the  image  to  the  closest,  next-smallest  power  of  2  dimention. 
If  n  =  5,  it  will  resize  to  2*5,  otherwise  it  will  use  2*4. 

This  is  necessary  for  creating  a  pyramid  to  represent  the  image  which  is 
at  least  4  levels  deep.  Bilinear  interpolation  is  used.  */ 
image  *  resize (image  *il ,  int  n) 

{ 

image  *newpic; 

int  xnew,  ynew,  preMc,  preMr,  index,  N,  x,  y; 

DPoint  color; 

if (n  ==  5)  N  =  32;  else  N  =  16; 


preMc  =  ( int) floor ( (double) il->xsize/N) ; 
preMr  =  (int) floor ( (double) il->ysize/N) ; 
xnew  =  preMc *N  +  1; 
ynew  =  preMr *N  +  1; 

newpic  =  (image  *)malloc (sizeof (image) ) ; 

newpic->pixel  =  (Pixel  *)malloc (xnew*ynew* sizeof (Pixel) ) ; 
newpic->xsize  =  xnew; 
newpic->ysize  =  ynew; 

for (y  =  0;  y  <  ynew;  y++)  { 

for (x  =  0;  x  <  xnew;  x++)  { 

index  =  y*xnew  +x; 

interpolate_pixel_bi (il,  xnew,  ynew,  x,  y,  color); 
newpic->pixel [index] [0]  =  color[0]; 
newpic- >pixel [index] [1]  =  color [1] ; 
newpic->pixel [index] [2]  =  color[2]; 

} 

} 

return (newpic) ; 

} 

/*  This  resizes  the  image  to  the  closest,  next-smallest  power  of  2  dimention. 
If  n  =  5,  it  will  resize  to  2*5,  otherwise  it  will  use  2*4. 

This  is  necessary  for  creating  a  pyramid  to  represent  the  image  which  is 
at  least  4  levels  deep.  Nearest  Neighbor  interpolation  is  used.  */ 
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image  *  res i z e_nn { image  *il,  int  n) 

{ 

image  *newpic; 

int  xnew,  ynew,  preMc,  preMr,  index ,  N,  x,  y; 

DPoint  color; 

if (n  ==  5)  N  =  32;  else  N  =  16; 

preMc  =  (int) floor ( (double) il->xsize/N) ; 
preMr  =  (int) floor ( (double) il->ysize/N) ; 
xnew  =  preMc *N  +  1; 
ynew  =  preMr *N  +  1; 

newpic  =  (image  *)malloc (sizeof (image) ) ; 

newpic->pixel  =  (Pixel  *)malloc (xnew*ynew*sizeof (Pixel) ) ; 
newpic- >xsize  =  xnew; 
newpic->ysize  =  ynew; 

for(y  =  0;  y  <  ynew;  y++)  { 

for (x  =  0;  x  <  xnew;  x++)  { 

index  =  y*xnew  +  x; 

interpolate_pixel_nn (il,  xnew,  ynew,  x,  y,  color); 
newpic->pixel [index] [0 ]  =  color[0]; 
newpic->pixel [index] tl]  =  color [1] ; 
newpic->pixel [index] [2 ]  =  color[2]; 

} 

} 

return (newpic) ; 


/*  This  resizes  the  image  to  a  defined  x  and  y  size;  defined  by  xnew 
and  ynew.  Bilinear  interpolation  is  used.  */ 
void  resize_def ine ( image  *il,  int  xnew,  int  ynew) 

{ 

Pixel  *pic; 
int  index,  x,  y; 

DPoint  color; 

pic  =  (Pixel  *)malloc (xnew*ynew*sizeof (Pixel) ) ; 

for(y  =0;  y  <  ynew;  y++)  { 

for (x  =  0;  x  <  xnew;  x++)  { 

index  =  y*xnew  +  x; 

interpolate^ ixel_bi (il,  xnew,  ynew,  x,  y,  color); 
pic [index] [0]  =  color[0]; 
pic [index] [1]  =  color[l]; 
pic [index] [2 ]  =  color[2]; 

> 

> 

free (il->pixel) ; 
il->pixel  =  pic; 
il->xsize  =  xnew; 
il->ysize  =  ynew; 
return ; 

} 

/♦Nearest  neighbor  interpolation  * / 

void  interpolate_pixel_nn ( image  *  g,  int  xnew,  int  ynew,  int  x,  int  y,  DPoint  color) 

{ 

int  fx,  fy,  index; 
double  xval,  yval; 

xval  =  (double) (x*g->xsize) / (double) xnew; 
yval  =  (double) (y*g->ysize) / (double) ynew; 

fx  =  (int) round (xval) ; 
f y  =  ( int ) round (yval ) ; 
index  =  fy*g->xsize  +  fx; 


color [0]  =  g->pixel [index] [0] ; 


color [1]  =  g->pixel  [index]  [1]  ; 
color[2]  =  g->pixel [index]  [2]  ; 
return; 

> 

/*bilinear  interpolation  */ 

void  interpolate _jpixel_bi (image  *  g,  int  xnew,  int  ynew,  int  x,  int  y,  DPoint  color) 

{ 

double  il,  i2f  i3 ,  i4; 
int  fx,  fy,  cx,  cy; 

double  xval,  yval,  slf  s2,  s3 ,  s4,  error  =  .0000000001? 

xval.  =  (x*g->xsize) /xnew; 
yval  =  (y*g->ysize) /ynew? 

fx  =  (int) floor (xval +error) ; 
cx  =  (int) ceil (xval) ; 

if (fx  ==  cx  &&  cx  !=  (g->xsize  -1))  cx++? 
fy  =  (int) floor (yval+error) ; 
cy  =  ( int) ceil (yval) ; 

if (fy  ==  cy  &&  cy  !=  (g->ysize  -1))  cy++? 

si  =  (xval  -  fx)/2; 
s2  =  0.5  -  si; 
s3  =  (yval  -  fy)/2; 
s4  =  0.5  -  s3 ; 


11  =  g->pixel [fy*g->xsize  + 

12  =  g->pixel  [fy*g*”>xsize  + 

13  =  g->pixel [cy*g->xsize  + 

14  =  g->pixel [cy*g->xsize  + 
color [0 ]  =  (sl*il  +  s2*i2  + 

11  =  g->pixel [fy*g->xsize  + 

12  =  g->pixel [fy*g->xsize  + 

13  =  g->pixel [cy*g->xsize  + 

14  =  g->pixel [cy*g->xsize  + 
color[l]  =  (sl*il  +  s2*i2  + 


fx]  [0]  ; 
cx] [0] ; 
fx] [01 ; 
cx] [0] ? 
s3*i3  + 
fx] [1]; 
cx] [1] ; 
fx] [1]; 
cx] [13; 
s3*i3  + 


s4*i4) ; 


s4*i4) ? 


} 


11  =  g->pixel [fy*g~>xsize  + 

12  =  g->pixel [fy*g->xsize  + 

13  =  g->pixel [cy*g->xsize  + 

14  =  g->pixel [cy*g->xsize  + 
color [2 ]  =  (sl*il  +  s2*i2  + 
return; 


fx]  [2]  ; 
cx] [2] ; 
fx]  [2] ; 
cx] [2] ; 
s3*i3  + 


s4*i4) ; 


#endif 


